小麥加工鏈中重金屬鎘含量的深度網(wǎng)絡(luò)預(yù)測模型

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

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

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

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

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

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

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

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

1 材料與方法

1.1 材料

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

1.2 設(shè)備

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

1.3 方法

1.3.1 小麥加工鏈及關(guān)鍵節(jié)點(diǎn)

小麥加工鏈?zhǔn)侵笍脑湹叫←湻鄣乃猩a(chǎn)環(huán)節(jié)，簡化示意圖如圖1所示。該加工鏈共有清理、潤麥、皮磨等10 個(gè)加工環(huán)節(jié)。原糧委托某公司按圖1進(jìn)行加工，并對(duì)各加工環(huán)節(jié)進(jìn)行抽檢。抽檢方法為：每隔2 h對(duì)全加工鏈抽檢1 次，共抽取16 次，每次在各環(huán)節(jié)分別抽取3～4 份（每份1 kg）平行樣品。為便于表示各環(huán)節(jié)的樣本，設(shè)定原糧、清理、潤麥、皮磨一（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部分：電熱原子吸收光譜法測定鎘含量》測定鎘含量（以鎘元素計(jì)），單位為mg/kg。

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

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

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

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

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

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

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

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

圖2 GRU模型結(jié)構(gòu)圖Fig. 2 Structure diagram of GRU model

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1.4.4 預(yù)測效果評(píng)價(jià)標(biāo)準(zhǔn)與模型對(duì)比驗(yàn)證

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

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

2 結(jié)果與分析

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

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

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

2.2 不同模型預(yù)測的性能對(duì)比

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

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

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

2.3 模型預(yù)測結(jié)果

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

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

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

3 結(jié) 論

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

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