基于BP神經網絡的污染場地土壤重金屬和PAHs含量預測

任加國， 龔 克,， 馬福俊， 谷慶寶， 武倩倩

1.山東科技大學地球科學與工程學院, 山東 青島 266590 2.中國環境科學研究院, 環境基準與風險評估國家重點實驗室, 北京 100012

工業活動產生的廢氣、廢液和廢渣排放及泄漏[1]導致大量重金屬和多環芳烴(polycyclic aromatic hydrocarbons，PAHs)進入土壤環境，造成土壤污染，其中以焦化廠、煉鋼廠、金屬加工廠及其周邊地區土壤污染尤為顯著[2-3]. 重金屬和PAHs具有性質穩定、難降解和毒性強的特點，容易被土壤吸附并不斷累積[4]，進而威脅人類健康[5-6]. 從保障人類生產生活安全角度出發，實現土壤環境質量評價和安全管理需要對土壤中重金屬和PAHs含量進行系統、全面地檢測. 然而，受土壤檢測成本和項目周期等因素制約，污染場地土壤經常存在檢測數據缺失的現象，如何利用有限的檢測數據獲得更全面的信息成為當前研究熱點.

人工神經網絡是模仿人腦結構和功能的運算模型，具有聯想記憶、分類識別、優化計算以及非線性映射等近似人腦的基本功能[7]. BP神經網絡是一種按照誤差逆傳播的人工神經網絡[8]，能夠處理已知條件與結果之間無明確關系的數據，通過在條件與結果之間建立一定的映射關系，而不需要在構建網絡之前確定映射的數學方程[9]，從不完整的樣本中提取信息特征對問題進行預測評估[10]. 目前，BP神經網絡對土壤污染物含量預測主要是針對空間尺度大區域的表層土壤[11-13]，而對小區域土壤污染物垂向空間分布預測的相關研究較少. 與大區域相比，小區域土壤污染受到其他因素(自然、生產因素)的影響較小，但土壤自身理化性質及其與污染物間的相互作用對污染物空間分布具有顯著影響，且土壤污染物空間分布的檢測能力及成本要求較高. 因此，通過BP神經網絡實現小區域場地土壤污染物的空間分布預測具有重要的現實意義.

某金屬加工廠是新中國成立后的一個重點綜合性有色金屬加工、科研、生產基地. 由于該廠對環保工作重視度不夠，設施落后，生產過程中產生大量的粉塵、酸洗廢液等，導致地塊內土壤受到重金屬和PAHs不同程度的污染. 以該場地為研究對象，對場地不同點位及深度的土壤樣品進行重金屬和PAHs含量檢測，分析不同污染物之間的關聯性，建立以監測點空間參數和已知土壤污染物含量數據為輸入值、土壤重金屬和PAHs含量數據為輸出值的映射關系，通過構建BP神經網絡模型，預測缺失的重金屬和PAHs含量數據，并對BP神經網絡預測效果進行評價，以期為土壤污染的空間分析和評價方法提供參考.

1 材料與方法

1.1 研究區概況及采樣點布置

研究區為某有色金屬加工廠生產區域，占地面積約6.5×105m2，自20世紀50年代開始，廠區一直從事有色金屬加工生產，主要產品包括銅、鎳、鈦合金和其他冶金制品等. 常年生產活動對地塊土壤造成了一定污染. 場地金屬加工車間(熔鑄、板帶、管棒車間等)主要沿南北方向分布在廠區中部，輔助系統和供暖車間分布在廠區西側. 場地地層自上而下主要為雜填土、粉質黏土、砂土和礫砂土，含水層埋藏較深. 按照HJ 25.2—2019《建設用地土壤污染風險管控和修復監測技術導則》相關要求，采用系統網格布點法，不超過 1 600 m2布設1個采樣點，對于存在污染區域進行加密布點，不超過400 m2布設1個采樣點. 場地總共設置379個采樣點，分別采集0.5～11.0 m不同深度的土壤樣品，其中表層土壤采樣深度為0～0.5 m，其他層位每1 m采集1個土壤樣品，共采集 1 832 個土壤樣品. 采樣點布設見圖1.

圖1 采樣點位置示意Fig.1 The location of sampling sites

1.2 土壤樣品測定

采集的土壤樣品剔除樹枝石塊等雜質，經冷凍干燥、混勻風干后研磨破碎過100目(0.15 mm)鋼篩后保存. 土壤pH和含水率分別利用pH計和含水率測定儀進行測定[14-15]；采用HCl-HNO3-HF微波密閉消解技術進行土壤樣品消解，使用電感耦合等離子體原子發射光譜法(ICAP RQ，Thermo Fisher Scientific，美國)測定土壤中w(Zn)、w(Cu)、w(Ni)、w(Cd)、w(Hg)、w(Cr)[16]，使用原子熒光法(AFS 2100，北京海光儀器有限公司)測定土壤中w(As)[17].

為保證分析樣品的準確性，采用現場平行樣、空白樣品、實驗室平行樣和加標樣品進行質量控制，國家土壤標準物質測試結果均在土壤成分的認定值范圍內，重金屬的加標回收率范圍為81%～130%，相對偏差控制在0～12%之間；PAHs的加標回收率范圍為55%～104%，相對偏差范圍控制在1%～14%之間. 測試結果均符合相應標準方法質量控制與保證要求.

1.3 土壤污染物關聯性分析

土壤污染物含量是多種因素共同影響的結果，其作用機理較為復雜，利用神經網絡進行污染物含量預測之前，需要考慮不同污染物之間的交互作用[19]，對模型的輸入參數進行選取，因此該文通過多元統計方法〔相關性分析、主成分分析(principal component analysis, PCA)和聚類分析〕分析土壤中重金屬和PAHs污染物間的關聯性，為神經網絡輸入參數的選取提供依據.

關聯性分析中重金屬污染物輸入參數為w(As)、w(Cr)、w(Zn)、w(Cu)、w(Pb)、w(Ni)、w(Cd)，PAHs污染物輸入參數為w(BaP)、w(DBA)、w(BkF)、w(BbF)、w(BaA)、w(Nap)、w(Chr). 相關性分析選擇雙變量Pearson法，顯著性檢驗選擇雙尾檢驗[20]；主成分分析選擇KMO和Bartlett法進行球形度檢驗，利用最大方差法獲得旋轉因子，旋轉方法為Kaiser標準化的正交旋轉法，旋轉在迭代5次后收斂[21-22]；聚類分析選擇按照組間連接系統聚類，選擇Pearson相關性作為度量標準.

關聯性顯著的元素可能具有同源性，能夠提高預測精度[23-24]，但是即使部分污染物關聯性較弱，二者之間可能存在非線性相關性，仍可能對預測精度有積極影響. 因此，該文中關聯性較弱的污染物數據仍作為神經網絡預測模型的輸入因子進行訓練，并設置將關聯性較弱因子去除后的輸入因子預測模型作為對照樣本，探究關聯性弱的因子對BP神經網絡模型預測精度的影響.

1.4 BP神經網絡模型構建及污染物含量預測

BP神經網絡是一種典型的多層前饋網絡[25]，其函數逼近、回歸計算等能力已被廣泛應用于環境科學研究領域，在預測土壤污染物的空間分布方面均取得了較好的效果[26-28].

該研究所選場地范圍較大，分析樣品數量多，受成本限制，個別采樣點的某些污染物指標未檢測，導致數據缺失. 為保證檢測數據的完整性，通過構建BP神經網絡對缺失數據進行預測：將樣本中不含缺失數據的因子作為模型的輸入條件，含有缺失數據的因子作為輸出條件，利用樣本中的已知數據訓練神經網絡，在訓練達到要求后，將缺失樣本的已知數據輸入模型，輸出值即為缺失數據的預測值[29]. 研究區共獲得1 691組土壤樣本數據，隨機抽取30組作為驗證樣本，在其余樣本中隨機選取50、200、800和1 661組(剩余全部樣本)作為訓練樣本，其中，w(Cu)、w(Pb)和w(DBA)、w(BkF)、w(BbF)、w(BaA)、w(Nap)、w(Chr)作為缺失數據，通過構建的神經網絡對缺失數據進行預測，并與驗證樣本的實測數據進行對比，評價模型的預測精度.

重金屬含量預測選取采樣點空間坐標參數、土壤pH、含水率(ω)、w(As)、w(Cr)、w(Zn)、w(Ni)、w(Cd)作為輸入參數，w(Pb)和w(Cd)作為輸出參數；PAHs含量預測選取采樣點空間坐標參數、土壤pH、含水率(ω)和w(BaP)作為輸入參數，w(DBA)、w(BkF)、w(BbF)、w(BaA)、w(Nap)、w(Chr)作為輸出參數.

構建的BP神經網絡拓撲結構見圖2，神經網絡每個節點代表一種特定的輸出函數，稱為激活函數[30]，每兩個節點的連接代表權重，通過不斷學習來調整權重值[31]，該文采用的學習速率為0.05. 構建的雙層隱含層單元數分別為20和15個. 選取正切函數tansig作為隱含層的傳遞函數，線性函數purelin作為輸出層的傳遞函數，共軛梯度函數trainscg作為訓練函數對樣本進行訓練[32]. 重金屬和PAHs指標最大迭代次數分別設置為20 000和 10 000 次，通過反復迭代，最終確定權值和閾值，建立預測模型.

圖2 BP神經網絡拓撲結構Fig.2 Topological of BP neural network

BP神經網絡根據梯度下降法調節連接權值，使訓練誤差達到最小，為使輸入數據能夠適應傳遞函數，提高計算過程的收斂速度，需要對訓練數據進行歸一化處理[33]，處理后的變量取值范圍在[-1,1]之間，計算方法如式(1)所示：

(1)

式中，xi為輸入變量ti歸一化后的值，tmax和tmin分別為輸入變量X的最大值和最小值.

1.5 數據分析

所有試驗數據使用Microsoft Excel 2013進行處理，并采用SPSS 22軟件進行多元統計分析；利用Matlab 2014軟件編寫程序構建BP神經網絡的預測模型，采用OriginPro 9.0軟件進行模型精度分析.

2 結果與討論

2.1 土壤重金屬和PAHs含量特征

研究區土壤污染物含量統計結果見表1[34]. 結果顯示，除w(Chr)外，土壤中重金屬和PAHs均存在不同程度的超標現象. 土壤重金屬污染物中，除w(Cr)外，其他6種重金屬平均含量均高于當地土壤背景值，與GB 36600—2018《土壤環境質量 建設用地土壤污染風險管控標準》[35]中的風險篩選值對比，7種重金屬含量超標率表現為w(Ni)>w(Cu)>w(As)>w(Pb)>w(Zn)=w(Cd)>w(Cr)，最大超標倍數分別為43.93、8.45、5.80、8.05、5.03、17.70和2.98倍；土壤PAHs污染物中，w(BaP)的超標情況較為嚴重，局部超標倍數較高，與GB 36600—2018中土壤PAHs風險篩選值對比，除w(Chr)未超標外，其他6種PAHs含量超標率表現為w(BaP)>w(Chr)>w(DBA)>w(BbF)=w(BaA)>w(Nap)>w(BkF)，最大超標倍數分別為120.00、28.00、30.55、30.36、3.02和1.22倍.

變異系數(CV)可以反映土壤污染物的變異程度，該值越大，表明該污染物受人為活動的影響越強[36]. 研究區土壤重金屬和PAHs的變異系數見表1[34]，其中重金屬Zn、Cu、Pb、Ni和Cd的變異系數范圍為2.75～8.12，屬于強變異(CV≥1.0)；As和Cr的變異系數分別為0.43和0.82，屬于中等變異(0.1≤CV<1.0)；PAHs變異系數為3.49～6.84，屬于強變異[37]. 由此可見，場地重金屬和PAHs污染物受到較強的人為污染源影響.

表1 研究區土壤重金屬和PAHs的含量統計值[34]

2.2 土壤污染物關聯性分析

2.2.1重金屬污染物

相關性分析通常作為污染物同源鑒別的依據[38-39]. 從表2可以得出，Zn與Cu、Pb、Cd均具有較強相關性(R為0.579～0.673，P<0.01)，Cu與Ni相關性較好(R=0.519)，Cr與其他6種重金屬相關性均較差，其他重金屬間彼此相關性較弱，表明Cr與其他重金屬的來源均不同. 由于研究場地不同類型的生產車間較多，因此不同區域土壤重金屬污染來源差異性較大.

表2 研究區土壤7種重金屬含量的相關性

KMO檢驗值(0.607)和Bartlett球形度(ξ=4 045，P<0.05)分析結果表明，該研究中重金屬污染物含量數據適用主成分分析. 7種重金屬經PCA解析后〔見圖3(a)〕，獲得PC1、PC2、PC3和PC4共4個主成分，累計方差貢獻率為84.01%，各主成分方差貢獻率分別為23.54%、23.08%、21.82%和15.57%. 從因子載荷(F)來看，PC1以As(F=0.876)、Pb(F=0.735)和Zn(F=0.510)為主導，PC2以Cd(F=0.930)和Zn(F=0.734)為主導，PC3以Ni(F=0.935)和Cu(F=0.747)為主導，PC4以Cr(F=0.97)為主導. 值得注意的是，前2個主成分中，Zn的正載荷均較高，表明Zn污染來源較廣，且Zn與As、Pb、Cd具有同源性；Cr與其他6種重金屬來源均不同.

聚類分析將7種重金屬分為三簇〔見圖3(b)〕，第一簇為Zn、Pb、As和Cd；第二簇為Cu和Ni；第三簇為單獨Cr，簇距離越低，表示因子間關聯越顯著[40-41]. 由圖3(b)可見，第一簇中Zn與Pb、As、Cd的簇距離較小，表明Zn與這3種重金屬具有相似污染源；Cr與其他簇的距離均較大，表明Cr與其他6種重金屬的來源差異性較大.

圖3 研究區土壤中7種重金屬主成分和聚類分析Fig.3 Principal component diagram and dendrogram of seven heavy metals contents in soils of the study area

上述3種多元統計分析方法所得結果相近，綜合分析結果并結合場地生產歷史可以得出：不同生產車間的特征污染物不同，Zn、Pb、As和Cd是熔鑄車間的主要產品或原料，關聯性較好；Cu和Ni主要是管棒車間的主要產物，關聯性較好，這些重金屬污染主要源自不同生產線的原材料和產品；Cr僅在盤管車間涉及污染，且w(Cr)(平均值為54.9 mg/kg)與當地土壤背景值(57.7 mg/kg)接近，相關研究表明，Cr污染主要源于成土母質[42-43]，研究區Cr判斷為受少量工業影響的自然源，因此與其他6種重金屬關聯性較弱.

2.2.2PAHs污染物

研究區土壤中7種PAHs相關性分析結果見表3，除Nap外，BaP、DBA、BkF、BbF、BaA和Chr彼此間均呈顯著相關(P<0.01)，表明這6種PAHs具有相同或相似的來源；Nap與其他PAHs相關性較弱.

表3 研究區土壤中7種PAHs的相關性

KMO檢驗值(0.854)和Bartlett球形度(ξ=23 761，P<0.05)分析表明，PAHs含量數據適用主成分分析. 通過分析獲得2個主成分〔見圖4(a)〕，累計方差貢獻率為90.53%，PC1、PC2的方差貢獻率分別為71.73%和18.80%. 從因子載荷(F)來看，PC1以BaP(F=0.762)、DBA(F=0.874)、BkF(F=0.953)、BbF(F=0.968)、BaA(F=0.946)和Chr(F=0.950)為主導，PC2以Nap(F=0.973)為主導，這與PAHs相關性分析結果一致.

聚類分析也將7種PAHs分為兩簇〔見圖4(b)〕，第一簇為BaP、DBA、BkF、BbF、BaA、Chr，簇距離均小于10；第二簇為單獨Nap，進一步表明除Nap外的其他6種PAHs具有相同或相似污染源.

圖4 研究區土壤中7種PAHs主成分和聚類分析Fig.4 Principal component diagram and dendrogram of seven PAHs in soils of the study area

結合廠區生產歷史并綜合3種分析方法結果表明，除Nap外，BaP、DBA、BkF、BbF、BaA和Chr彼此關聯性較好，其中BaP、BaA和Chr是燃煤型排放的標志污染物，BkF、BbF和DBA主要來自汽油或柴油的燃燒[44]，因此推斷這6種PAHs主要來自燃煤和交通混合源[45]；而Nap的產生可能與生物質不完全燃燒有關[46]，因此與其他6種PAHs關聯性較弱.

2.3 BP神經網絡訓練及檢驗

2.3.1不同訓練樣本數對神經網絡精度影響

按照1.4節建立的BP神經網絡模型進行訓練，得到不同訓練樣本數量對訓練精度的影響如圖5所示. 由圖5可以看出，隨著訓練樣本數量的增加，神經網絡訓練誤差(mean squared error，MSE)逐漸減小，并且訓練結束時MSE更接近設置的目標誤差，表明隨著樣本數量的增加，訓練精度逐步提高；而當訓練樣本數據減至50個時，重金屬和PAHs相應BP模型在未達到迭代次數(20 000 和 10 000 次)時提前終止訓練，表明隨著迭代次數的增加，訓練誤差不再降低，訓練精度無法提高，停止訓練. 因此，適當增加訓練樣本數量能夠有效提高模型預測精度.

圖5 不同訓練樣本數量對BP神經網絡精度的影響Fig.5 The effect of different training sample size on the accuracy of BP neural network

2.3.2神經網絡訓練精度和誤差分析

選擇全部1 661個訓練樣本建立BP神經網絡模型，結果如圖6所示. 重金屬及PAHs樣本在訓練次數分別為 20 000 和 10 000 次時達到最佳效果，訓練樣本MSE分別為1.9×10-4和4.3×10-5，表明訓練收斂效果較好；同時，訓練樣本輸出值和目標值之間的相關系數(R)分別達到 0.990 11 和 0.995 19，R接近1，表明重金屬和PAHs樣本神經網絡訓練效果較好.

圖6 BP神經網絡訓練誤差曲線和線性回歸結果Fig.6 The error curve diagram and linear regression diagram of BP neural network training

利用訓練完成后的BP神經網絡模型對30組驗證樣本的缺失數據〔w(Cu)、w(Pb)和w(DBA)、w(BkF)、w(BbF)、w(BaA)、w(Nap)、w(Chr)〕進行預測，結果如圖7所示. 缺失數據預測結果與實測結果的R范圍為0.901～0.996，預測曲線與實測曲線重合度較好. 將BP神經網絡模型預測結果與實測結果進行對比分析，結果(見表4)顯示，驗證樣本的決定系數(R2)范圍為0.812～0.993，表明該模型預測擬合度效果較高. 模擬效率系數(Nash-Sutcliffe efficiency coefficient，NSE)通常是衡量模擬結果可靠性的重要指標[29]，取值范圍為(-∞，1]，當NSE為正值時，表明模擬結果可信，且該值越接近1，表示模型匹配程度越好. 由表4可見：該研究中NSE范圍為0.779～0.959，除w(DBA)外，其他驗證樣本的NSE均大于0.8，模擬效果較好，這與其他相關研究結果[47]相似；同時，驗證樣本的均方根誤差(RMSE)和平均絕對誤差(MAE)均較小. 綜上，BP神經網絡模型可靠性高，能夠較好地預測研究區土壤重金屬和PAHs含量.

圖7 BP神經網絡預測值與實測值結果散點圖Fig.7 The scatter diagram of predicted values and measured values of BP neural network

表4 BP神經網絡預測精度和誤差分析

為了探究污染物關聯性對BP神經網絡模型預測精度的影響，根據污染物的關聯性分析結果，該研究去除了關聯性較弱的因子〔w(Cd)、w(Cr)〕，選取重金屬中關聯性較強的因子〔(w(As)、w(Cr)、w(Zn)、w(Ni)〕和采樣點空間參數以及理化性質作為輸入變量，w(Cu)、w(Pb)作為輸出變量，在其他參數不變的條件下構建對照模型，該對照模型的預測結果與選取全部因子作為輸入變量構建的模型預測結果對比見表5. 由表5可以看出，與去除關聯性較弱的因子相比，選取全部因子作為輸入變量建立的預測模型整體擬合度更好，其精度分析指標(R2)相比于對照模型提高了0.129，NSE提高了0.134，誤差分析指標RMSE和MAE分別降低了0.036和0.031. 由此可見，相比于僅使用關聯性強的因子作為輸入參數，考慮不同影響因子構建的BP神經網絡模型預測更加準確可靠. 推測原因可能是，關聯性較弱因子之間仍有非線性相關性，但這種關聯難以通過具體數學方程計算，而神經網絡不需要確定具體的映射方程，因此將這些因子作為輸入變量構建神經網絡，能夠和輸出端建立較好的映射關系，從而進一步提高預測模型的精度.

表5 不同輸入變量的預測模型精度對比

3 結論

a) 通過分析某金屬加工廠生產區內土壤污染特性，結果表明，除PAHs的Chr外，重金屬和PAHs均存在不同程度的超標現象，其中重金屬Ni、As、Cu及PAHs中BaP超標較為嚴重；重金屬和PAHs變異系數較高，該區域土壤受到較強的人為污染源影響.

b) 多元統計分析結果顯示，污染土壤中重金屬Zn與Pb、As、Cd關聯性均較好，Cu與Ni關聯性較好，而Cr與其他6種重金屬關聯性較弱，Zn、Pb、As、Cd、Cu和Ni污染主要源于不同生產線的原材料和產品，Cr為受少量工業影響的自然源；PAHs中除Nap外，BaP、DBA、BkF、BbF、BaA和Chr彼此關聯性較好，6種關聯性好的PAHs主要來自燃煤和交通混合源，Nap則源于生物質不完全燃燒.

c) 構建的BP神經網絡模型訓練效果較好；驗證樣本的預測值與實測值之間的誤差分析顯示，各污染物含量的決定系數(R2)范圍為0.812～0.993，模擬效率系數(NSE)范圍為0.779～0.959，均方根誤差(RMSE)和平均絕對誤差(MAE)均較小；構建的BP神經網絡模型準確可靠，能夠較好地預測研究區內土壤污染物含量，且關聯性較弱因子的輸入能進一步提高預測模型的精度.