一種分組模式下的土壤重金屬含量預(yù)測模型

呂鑫濤，張 聰，曹文琪

（武漢輕工大學(xué)數(shù)學(xué)與計算機(jī)學(xué)院，湖北武漢 430023）

0 引言

土壤是人類生活、社會發(fā)展和自然進(jìn)化的重要資源。要改善地球生態(tài)環(huán)境和人類生活環(huán)境，對土壤進(jìn)行整治必不可少，但土壤中的重金屬是非常棘手的一種污染物，難以被微生物降解，影響農(nóng)作物生長，導(dǎo)致農(nóng)作物產(chǎn)量下降；通過人體接觸或食用最終進(jìn)入人體，影響健康［1-5］。因此，對土壤重金屬進(jìn)行研究成為環(huán)境污染研究的重要內(nèi)容之一。

土壤重金屬含量預(yù)測研究提出了多種預(yù)測模型，對土壤屬性作了相關(guān)分析，如陳飛香等［6］以廣東增城市為例，采用隨機(jī)采樣方法對土壤中Cr 的含量進(jìn)行分析，通過RBF 網(wǎng)絡(luò)對Cr 含量插值進(jìn)行預(yù)測誤差分析，證實RBF 網(wǎng)絡(luò)插值方法克服了克里格插值的平滑效應(yīng)，有更好的預(yù)測效果；Ser?geev 等［7］提出一種人工神經(jīng)網(wǎng)絡(luò)（ANN），模擬非線性大規(guī)模趨勢，用統(tǒng)計方法對殘差建模預(yù)測土壤重金屬，證實神經(jīng)網(wǎng)絡(luò)預(yù)測重金屬空間分布的可行性；張鈺等［8］將新疆克拉瑪依人工碳匯林作為采樣區(qū)域，通過BP 神經(jīng)網(wǎng)絡(luò)和GIS技術(shù)建模，對土壤中Cu、Zn、Fe 等5 種金屬含量的分布進(jìn)行分析，證實雙隱層BP 神經(jīng)網(wǎng)絡(luò)能夠滿足預(yù)測精度并能分析土壤中重金屬含量的空間分布；楊勇等［9］用時空克里格模型對武漢市青山區(qū)土壤中重金屬含量在時間和空間進(jìn)行了建模預(yù)測，在土壤污染防控治理問題上具有參考作用；呂杰等［10］以陜西金堆城礦區(qū)尾礦庫為例，通過實驗室化學(xué)分析和小波神經(jīng)網(wǎng)絡(luò)對土壤中Cu 的含量進(jìn)行了分析，利用ASD 光譜儀的遙感數(shù)據(jù)對土壤重金屬做出評價，證實小波神經(jīng)網(wǎng)絡(luò)反演土壤Cu 元素含量的有效性；Cao 等［11-12］以小波神經(jīng)網(wǎng)絡(luò)和徑向基神經(jīng)網(wǎng)絡(luò)為主分別提出了協(xié)同復(fù)合神經(jīng)網(wǎng)絡(luò)模型和深度復(fù)合模型，通過對神經(jīng)網(wǎng)絡(luò)的初始權(quán)值與閾值進(jìn)行優(yōu)化，并且對土壤重金屬含量進(jìn)行預(yù)測分析，提升了模型的預(yù)測精度；Eid 等［13］運(yùn)用污水污泥（sew?age sludge）改良過的土壤種植作物，可以提高模型預(yù)測的準(zhǔn)確性。黃趙麟等［14］構(gòu)建了源匯模型（BP-S）、空間分異模型（BP-K）和改進(jìn)的多因素綜合模型（BP-SK）模擬預(yù)測了5 種重金屬Cd、Pb、Cr、Cu 和Zn 含量，證實了BP-SK 模型在人為干擾大的區(qū)域中重金屬空間分布預(yù)測比BP-K 模型更實用；Zahida 等［15］對卡拉奇城市土壤進(jìn)行污染評價，利用累積頻率分布（CDF）曲線估算了重金屬地球化學(xué)基準(zhǔn)濃度；Pandit 等［16］通過偏最小二乘法（PLSR）建立模型反映土壤重金屬含量與高光譜反射率之間的關(guān)系，證實光譜反射率在反演Pd 等重金屬濃度時有著較高的預(yù)測精度，也說明光譜反射技術(shù)能夠繪制重金屬的空間分布；Vinod 等［17］用了回歸建模，主成分分析（PCA）和積累營養(yǎng)成分（ANE）方法對綜合工業(yè)廢水灌溉土壤重金屬進(jìn)行分析，提出了以土壤pH 值和重金屬含量作為自變量的模型。上述對土壤重金屬環(huán)境污染進(jìn)行評估，以及對土壤時間和空間特征分布的研究都取得了不錯的預(yù)測效果。本文提出一種基于分組教學(xué)優(yōu)化算法［18］（Group Teaching Optimization Algorithm，GTOA）和BP 神經(jīng)網(wǎng)絡(luò)預(yù)測模型，能夠有效避免BP 神經(jīng)網(wǎng)絡(luò)在訓(xùn)練時收斂過慢問題，提升預(yù)測精度。

1 基礎(chǔ)理論

1.1 分組教學(xué)優(yōu)化算法

分組教學(xué)優(yōu)化算法相較于傳統(tǒng)的智能優(yōu)化算法，最顯著的特點(diǎn)是需要調(diào)整的參數(shù)很少，只有種群大小以及迭代次數(shù)兩個參數(shù)，這種少量參數(shù)設(shè)置能極大發(fā)揮算法的優(yōu)化性能。算法流程主要分為能力分組階段、教師教學(xué)階段、學(xué)生學(xué)習(xí)階段、教師分配階段。

（1）能力分組階段。對于一個班級的知識分布，不失一般性可以假設(shè)其服從正態(tài)分布。教學(xué)不僅要提高班級的平均水平，還要考慮怎么減少標(biāo)準(zhǔn)差，使整個班級的知識分布更加合理且優(yōu)秀。要實現(xiàn)這一教學(xué)目標(biāo)，教師就需要針對不同的學(xué)生制定相應(yīng)的教學(xué)計劃。將學(xué)生分成兩個小組，接受知識能力強(qiáng)的稱為優(yōu)秀學(xué)生，接受知識能力弱的稱為一般學(xué)生。

（2）教師教學(xué)階段。對于出色學(xué)生組，通過下面方式獲取知識：

其中，t是當(dāng)前的迭代次數(shù)，是i學(xué)生在t時刻時的知識，Tt是時刻t時教師所擁有的知識，Mt是在t時刻小組的平均水平知識，F(xiàn)是決定教師教學(xué)成果的教學(xué)因素，F(xiàn)值可以是1 或2，a,b,c都是［0，1］之間的隨機(jī)數(shù)。

對于一般學(xué)生，獲取知識的方式如下：

其中，d是［0，1］之間的隨機(jī)數(shù)。

（3）學(xué)生學(xué)習(xí)階段。學(xué)生學(xué)習(xí)階段又可分成兩個階段，除開教師教學(xué)時間，學(xué)生可以在課余時間通過自學(xué)以及與其他優(yōu)秀同學(xué)交流來獲取知識，這一過程用下式來概括：

其中,e和g是［0，1］之間的隨機(jī)數(shù)表示在t時刻時學(xué)生j在學(xué)習(xí)階段學(xué)到的知識，學(xué)生j∈{1,2,…,i-1,i+1,…,N}是隨機(jī)的。向其他學(xué)生交流學(xué)習(xí)和自學(xué)階段學(xué)習(xí)分別是式（4）中的第2 項和第3 項。

（4）教師分配階段。為保證算法的收斂速度，優(yōu)秀學(xué)生組和一般學(xué)生組共用一位教師，而對教師的選擇則受到灰狼優(yōu)化算法（Grey Wolf Optimizer，GWO）中保留3 個最優(yōu)的思想影響。具體的教師分配可根據(jù)下式選擇：

GTOA 算法步驟如下：①初始化相關(guān)參數(shù)，隨機(jī)產(chǎn)生含有N 個個體的學(xué)生種群；②計算個體的適應(yīng)度值，選出最優(yōu)解，并且更新函數(shù)迭代次數(shù)；③判斷是否達(dá)到最大迭代次數(shù)，達(dá)到了即終止循環(huán)并輸出最優(yōu)解，否則跳轉(zhuǎn)至第④步；④拿到3 個最優(yōu)解，根據(jù)公式（5）選出教師；⑤基于個體的適應(yīng)度值大小分組，將適應(yīng)度值大的一半個體組成優(yōu)秀學(xué)生組，將適應(yīng)度值小的另一半個體組成一般學(xué)生組，兩組學(xué)生共用一位教師；⑥對于優(yōu)秀學(xué)生組根據(jù)式（1）、式（2）和式（4）更新種群；對于一般學(xué)生組，根據(jù)式（3）和式（4）分別執(zhí)行教師教學(xué)階段和學(xué)生學(xué)習(xí)階段；⑦將兩個新的種群合并組成一個新的種群，計算種群的適應(yīng)度值，選取最優(yōu)解，并且更新迭代次數(shù)，然后返回到第③步。

1.2 BP 神經(jīng)網(wǎng)絡(luò)

BP 神經(jīng)網(wǎng)絡(luò)是在人工神經(jīng)網(wǎng)絡(luò)基礎(chǔ)上加入反向傳播算法的多層前饋神經(jīng)網(wǎng)絡(luò)。BP 神經(jīng)網(wǎng)絡(luò)特點(diǎn)是信號正向傳播，通過隱藏層和輸出層的計算得到誤差值，然后將誤差一層一層地反向傳播，更新網(wǎng)絡(luò)中的權(quán)值和閾值，達(dá)到最優(yōu)的期望值［19-20］。BP 神經(jīng)網(wǎng)絡(luò)（以一層隱含層為例）結(jié)構(gòu)如圖1 所示。

Fig.1 Three layer BP neural network structure圖1 三層BP 神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)

X1,…,Xn是網(wǎng)絡(luò)的輸入?yún)?shù)，Y1,…,Yn是網(wǎng)絡(luò)的輸出參數(shù)。初始化輸入層個數(shù)n，隱藏層個數(shù)l，輸出層個數(shù)m，輸入?yún)?shù)經(jīng)過輸入層到達(dá)隱藏層公式如下：

其中，f是隱藏層激活函數(shù)，可以使線性數(shù)據(jù)集變成非線性數(shù)據(jù)集。Wij和bj分別是輸入層到隱藏層之間的權(quán)值和閾值。由隱藏層到輸出層的計算公式如下：

其中,f是輸出層的激活函數(shù)，與隱藏層是同一個函數(shù)。Wjk和bk分別是隱藏層到輸出層之間的權(quán)值和閾值。Yk是網(wǎng)絡(luò)的預(yù)測輸出值。得到輸出值后計算誤差，公式如下：

將誤差反向傳播到隱藏層，對隱藏層和輸出層之間的權(quán)值和閾值進(jìn)行更新，公式如下：

為使誤差變小，再將誤差傳播至輸入層，如此反復(fù)迭代訓(xùn)練，判斷網(wǎng)絡(luò)是否達(dá)到最大迭代次數(shù)或者達(dá)到最小誤差條件。若沒有結(jié)束，則繼續(xù)使用更新后的權(quán)值與閾值進(jìn)入網(wǎng)絡(luò)循環(huán)，直至滿足網(wǎng)絡(luò)結(jié)束條件。

1.3 分組模式模型（GTOA-BP）

文獻(xiàn)［18］提出的GTOA 算法在約束工程設(shè)計優(yōu)化問題的廣義求解方法中有著很好的結(jié)果，比如焊接梁設(shè)計問題、管柱設(shè)計問題、壓力容器設(shè)計問題以及減速機(jī)設(shè)計問題，GTOA 優(yōu)秀的全局搜索能力能使BP 神經(jīng)網(wǎng)絡(luò)中的性能函數(shù)實現(xiàn)快速收斂，得到GTOA-BP 中的最優(yōu)權(quán)值與閾值，從而減少神經(jīng)網(wǎng)絡(luò)訓(xùn)練過程中產(chǎn)生的誤差，提高預(yù)測精度。GTOA 算法由學(xué)生個體組成的種群與網(wǎng)絡(luò)中的權(quán)值相對應(yīng)，適應(yīng)度函數(shù)是網(wǎng)絡(luò)訓(xùn)練常用的性能函數(shù)均方誤差（MSE）。模型具體步驟如下：

（1）對數(shù)據(jù)中農(nóng)作物按照類型進(jìn)行編號，然后將數(shù)據(jù)中樣本的經(jīng)度、緯度、海拔以及農(nóng)作物編號4 個特征值作為神經(jīng)網(wǎng)絡(luò)輸入?yún)?shù)，以重金屬As（mg/kg）含量作為輸出參數(shù)，對輸入?yún)?shù)和輸出參數(shù)最大、最小作歸一化處理。

（2）確定BP 神經(jīng)網(wǎng)絡(luò)中隱藏層個數(shù)、激活函數(shù)、學(xué)習(xí)率、性能函數(shù)、網(wǎng)絡(luò)訓(xùn)練次數(shù)以及GTOA 算法中種群大小和教學(xué)因素參數(shù)。

（3）初始化由權(quán)值與閾值個體組成的種群：

其中,N是生成種群大小，D是輸入層個數(shù)，即維度，u和l分別是權(quán)值和閾值的上下界，k是［0，1］之間的隨機(jī)數(shù)。

（4）計算種群適應(yīng)度值，按照適應(yīng)度值大小將種群分成優(yōu)秀種群和一般種群。

（5）從分組的種群中按式（5）的種群選取教師。

（6）更新種群位置。優(yōu)秀學(xué)生通過式（1）、式（2）、式（4），一般學(xué)生通過式（3）、式（4）完成教師教學(xué)階段和學(xué)生學(xué)習(xí)階段更新種群。

（7）重復(fù)步驟（4）-（6）直到滿足結(jié)束條件，得到最優(yōu)的權(quán)值和閾值。

具體流程如圖2 所示。

Fig.2 GTOA-BP model flow圖2 GTOA-BP 模型流程

2 實驗與結(jié)果分析

2.1 數(shù)據(jù)來源

本實驗數(shù)據(jù)來自武漢市農(nóng)科院環(huán)安所。依據(jù)《土壤環(huán)境監(jiān)測技術(shù)規(guī)范》（HJ/T166-2004）和《土壤環(huán)境質(zhì)量農(nóng)業(yè)農(nóng)用地土壤環(huán)境污染管控》（GB15618-2018）的要求，對武漢市6 個新城區(qū)的周邊土壤中8 種重金屬進(jìn)行檢測，包括砷（As）、銅（Cd）、鉻（Cr）、銅（Cu）、鎳（Ni）、鉛（Pb）、鋅（Zn）和汞（Hg），通過GPS 定位儀得到采樣點(diǎn)的經(jīng)度、緯度和海拔，記錄了1161 個采樣點(diǎn)的種植物類型。本研究的所有實驗均在處理器為2.4 GHz 的四核Intel Core i5 以及操作系統(tǒng)為Windows 的電腦上進(jìn)行，編程語言為Python，編程軟件為PyCharm 2019 專業(yè)版。將數(shù)據(jù)中的農(nóng)作物按照類型進(jìn)行編號，把經(jīng)度、緯度和高度的特征值作為模型輸入?yún)?shù)，將8 種元素中的As 含量（mg/kg）作為實驗?zāi)Ｐ洼敵鰠?shù)，采用數(shù)據(jù)集中的500 組數(shù)據(jù)作為訓(xùn)練集，30 組數(shù)據(jù)作為測試集。

對數(shù)據(jù)進(jìn)行預(yù)處理公式如下：

其中，X是原始數(shù)據(jù)，Xmax和Xmin分別是原始數(shù)據(jù)中的最大值和最小值。

2.2 實驗結(jié)果分析

為了檢驗預(yù)測模型可行性，將GTOA-BP 模型分別與BP 神經(jīng)網(wǎng)絡(luò)模型和徑向基神經(jīng)網(wǎng)絡(luò)模型（Radial Basis Function Neyral Network，RBFNN）進(jìn)行對比，實驗中BP 模型和RBF 模型迭代次數(shù)為600，激活函數(shù)采用tanh 函數(shù)，隱藏層節(jié)點(diǎn)數(shù)為8，學(xué)習(xí)率為0.02。GTOA-BP 預(yù)測模型種群大小為30，教學(xué)因素F為2，3 種模型訓(xùn)練下的As 含量真實值與預(yù)測值對比如圖3-圖5 所示。

Fig.3 Comparison of As content in BP model圖3 BP 模型As 含量對比

Fig.4 Comparison of As content in RBF model圖4 RBF 模型As 含量對比

Fig.5 Comparison of As content in GTOA-BP model圖5 GTOA-BP 模型As 含量對比

將預(yù)測得到的值與真實值作差值運(yùn)算，然后得到差值與對應(yīng)點(diǎn)的真實值比值大小，得到30 個測試點(diǎn)的分布，如表1 所示。

Table 1 Distribution of ratio of prediction points表1 預(yù)測點(diǎn)比值大小分布情況

平均絕對誤差（MAE）、均方誤差（MSE）、平均絕對百分比誤差（MAPE）和對稱平均絕對百分比誤差（SMAPE）是4種可以評判預(yù)測效果的誤差指標(biāo)公式。計算公式如下：

結(jié)合圖3、圖4、圖5 以及表1 和 表2 的數(shù)據(jù)可知，GTOA-BP 模型預(yù)測的數(shù)據(jù)相較于BP 模型和RBF 模型更加貼合真實數(shù)據(jù)，預(yù)測數(shù)據(jù)與真實數(shù)據(jù)差距較小。從表1 可以看出，GTOA-BP 模型得到的數(shù)值比其他兩種模型可信度高，在比值小于10%時，GTOA-BP 有比其他兩種模型一倍的預(yù)測點(diǎn)。由表2 可以看出，BPNN 模型和RBFNN 模型的4種誤差指標(biāo)數(shù)值均高于GTOA-BP 模型的誤差指標(biāo)，說明GTOA-BP 模型在3 種模型中預(yù)測效果最好。

Table 2 Comparison of prediction errors表2 預(yù)測誤差對比

3 結(jié)語

本文提出的GTOA-BP 預(yù)測模型通過GTOA 對BP 神經(jīng)網(wǎng)絡(luò)中的權(quán)值與閾值進(jìn)行優(yōu)化，相較于BPNN 和RBFNN 模型在土壤重金屬含量預(yù)測效果方面有著更好的表現(xiàn)。通過分組模型對土壤中殘留物進(jìn)行預(yù)測，以此為基礎(chǔ)對土壤重金屬開展風(fēng)險評估工作，為開展土地預(yù)防治理提供技術(shù)支撐。后續(xù)可以嘗試將GTOA 運(yùn)用到其他不同神經(jīng)網(wǎng)絡(luò)中，如日常生活中。但GTOA 算法還未運(yùn)用于實際問題，算法還需進(jìn)行優(yōu)化以加快算法收斂速度。