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

呂鑫濤，張 聰，曹文琪

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

0 引言

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

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

1 基礎(chǔ)理論

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

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

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

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

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

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

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

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

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

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

GTOA 算法步驟如下：①初始化相關(guān)參數(shù)，隨機(jī)產(chǎn)生含有N 個(gè)個(gè)體的學(xué)生種群；②計(jì)算個(gè)體的適應(yīng)度值，選出最優(yōu)解，并且更新函數(shù)迭代次數(shù)；③判斷是否達(dá)到最大迭代次數(shù)，達(dá)到了即終止循環(huán)并輸出最優(yōu)解，否則跳轉(zhuǎn)至第④步；④拿到3 個(gè)最優(yōu)解，根據(jù)公式（5）選出教師；⑤基于個(gè)體的適應(yīng)度值大小分組，將適應(yīng)度值大的一半個(gè)體組成優(yōu)秀學(xué)生組，將適應(yīng)度值小的另一半個(gè)體組成一般學(xué)生組，兩組學(xué)生共用一位教師；⑥對(duì)于優(yōu)秀學(xué)生組根據(jù)式（1）、式（2）和式（4）更新種群；對(duì)于一般學(xué)生組，根據(jù)式（3）和式（4）分別執(zhí)行教師教學(xué)階段和學(xué)生學(xué)習(xí)階段；⑦將兩個(gè)新的種群合并組成一個(gè)新的種群，計(jì)算種群的適應(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)是信號(hào)正向傳播，通過(guò)隱藏層和輸出層的計(jì)算得到誤差值，然后將誤差一層一層地反向傳播，更新網(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ù)。初始化輸入層個(gè)數(shù)n，隱藏層個(gè)數(shù)l，輸出層個(gè)數(shù)m，輸入?yún)?shù)經(jīng)過(guò)輸入層到達(dá)隱藏層公式如下：

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

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

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

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

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

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

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

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

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

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

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

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

（6）更新種群位置。優(yōu)秀學(xué)生通過(guò)式（1）、式（2）、式（4），一般學(xué)生通過(guò)式（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 實(shí)驗(yàn)與結(jié)果分析

2.1 數(shù)據(jù)來(lái)源

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

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

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

2.2 實(shí)驗(yàn)結(jié)果分析

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

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

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

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

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

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

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

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

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

3 結(jié)語(yǔ)

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