任加國, 龔 克,, 馬???, 谷慶寶, 武倩倩
1.山東科技大學地球科學與工程學院, 山東 青島 266590 2.中國環(huán)境科學研究院, 環(huán)境基準與風險評估國家重點實驗室, 北京 100012
工業(yè)活動產生的廢氣、廢液和廢渣排放及泄漏[1]導致大量重金屬和多環(huán)芳烴(polycyclic aromatic hydrocarbons,PAHs)進入土壤環(huán)境,造成土壤污染,其中以焦化廠、煉鋼廠、金屬加工廠及其周邊地區(qū)土壤污染尤為顯著[2-3]. 重金屬和PAHs具有性質穩(wěn)定、難降解和毒性強的特點,容易被土壤吸附并不斷累積[4],進而威脅人類健康[5-6]. 從保障人類生產生活安全角度出發(fā),實現(xiàn)土壤環(huán)境質量評價和安全管理需要對土壤中重金屬和PAHs含量進行系統(tǒng)、全面地檢測. 然而,受土壤檢測成本和項目周期等因素制約,污染場地土壤經(jīng)常存在檢測數(shù)據(jù)缺失的現(xiàn)象,如何利用有限的檢測數(shù)據(jù)獲得更全面的信息成為當前研究熱點.
人工神經(jīng)網(wǎng)絡是模仿人腦結構和功能的運算模型,具有聯(lián)想記憶、分類識別、優(yōu)化計算以及非線性映射等近似人腦的基本功能[7]. BP神經(jīng)網(wǎng)絡是一種按照誤差逆?zhèn)鞑サ娜斯ど窠?jīng)網(wǎng)絡[8],能夠處理已知條件與結果之間無明確關系的數(shù)據(jù),通過在條件與結果之間建立一定的映射關系,而不需要在構建網(wǎng)絡之前確定映射的數(shù)學方程[9],從不完整的樣本中提取信息特征對問題進行預測評估[10]. 目前,BP神經(jīng)網(wǎng)絡對土壤污染物含量預測主要是針對空間尺度大區(qū)域的表層土壤[11-13],而對小區(qū)域土壤污染物垂向空間分布預測的相關研究較少. 與大區(qū)域相比,小區(qū)域土壤污染受到其他因素(自然、生產因素)的影響較小,但土壤自身理化性質及其與污染物間的相互作用對污染物空間分布具有顯著影響,且土壤污染物空間分布的檢測能力及成本要求較高. 因此,通過BP神經(jīng)網(wǎng)絡實現(xiàn)小區(qū)域場地土壤污染物的空間分布預測具有重要的現(xiàn)實意義.
某金屬加工廠是新中國成立后的一個重點綜合性有色金屬加工、科研、生產基地. 由于該廠對環(huán)保工作重視度不夠,設施落后,生產過程中產生大量的粉塵、酸洗廢液等,導致地塊內土壤受到重金屬和PAHs不同程度的污染. 以該場地為研究對象,對場地不同點位及深度的土壤樣品進行重金屬和PAHs含量檢測,分析不同污染物之間的關聯(lián)性,建立以監(jiān)測點空間參數(shù)和已知土壤污染物含量數(shù)據(jù)為輸入值、土壤重金屬和PAHs含量數(shù)據(jù)為輸出值的映射關系,通過構建BP神經(jīng)網(wǎng)絡模型,預測缺失的重金屬和PAHs含量數(shù)據(jù),并對BP神經(jīng)網(wǎng)絡預測效果進行評價,以期為土壤污染的空間分析和評價方法提供參考.
研究區(qū)為某有色金屬加工廠生產區(qū)域,占地面積約6.5×105m2,自20世紀50年代開始,廠區(qū)一直從事有色金屬加工生產,主要產品包括銅、鎳、鈦合金和其他冶金制品等. 常年生產活動對地塊土壤造成了一定污染. 場地金屬加工車間(熔鑄、板帶、管棒車間等)主要沿南北方向分布在廠區(qū)中部,輔助系統(tǒng)和供暖車間分布在廠區(qū)西側. 場地地層自上而下主要為雜填土、粉質黏土、砂土和礫砂土,含水層埋藏較深. 按照HJ 25.2—2019《建設用地土壤污染風險管控和修復監(jiān)測技術導則》相關要求,采用系統(tǒng)網(wǎng)格布點法,不超過 1 600 m2布設1個采樣點,對于存在污染區(qū)域進行加密布點,不超過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
采集的土壤樣品剔除樹枝石塊等雜質,經(jīng)冷凍干燥、混勻風干后研磨破碎過100目(0.15 mm)鋼篩后保存. 土壤pH和含水率分別利用pH計和含水率測定儀進行測定[14-15];采用HCl-HNO3-HF微波密閉消解技術進行土壤樣品消解,使用電感耦合等離子體原子發(fā)射光譜法(ICAP RQ,Thermo Fisher Scientific,美國)測定土壤中w(Zn)、w(Cu)、w(Ni)、w(Cd)、w(Hg)、w(Cr)[16],使用原子熒光法(AFS 2100,北京海光儀器有限公司)測定土壤中w(As)[17].
為保證分析樣品的準確性,采用現(xiàn)場平行樣、空白樣品、實驗室平行樣和加標樣品進行質量控制,國家土壤標準物質測試結果均在土壤成分的認定值范圍內,重金屬的加標回收率范圍為81%~130%,相對偏差控制在0~12%之間;PAHs的加標回收率范圍為55%~104%,相對偏差范圍控制在1%~14%之間. 測試結果均符合相應標準方法質量控制與保證要求.
土壤污染物含量是多種因素共同影響的結果,其作用機理較為復雜,利用神經(jīng)網(wǎng)絡進行污染物含量預測之前,需要考慮不同污染物之間的交互作用[19],對模型的輸入?yún)?shù)進行選取,因此該文通過多元統(tǒng)計方法〔相關性分析、主成分分析(principal component analysis, PCA)和聚類分析〕分析土壤中重金屬和PAHs污染物間的關聯(lián)性,為神經(jīng)網(wǎng)絡輸入?yún)?shù)的選取提供依據(jù).
關聯(lián)性分析中重金屬污染物輸入?yún)?shù)為w(As)、w(Cr)、w(Zn)、w(Cu)、w(Pb)、w(Ni)、w(Cd),PAHs污染物輸入?yún)?shù)為w(BaP)、w(DBA)、w(BkF)、w(BbF)、w(BaA)、w(Nap)、w(Chr). 相關性分析選擇雙變量Pearson法,顯著性檢驗選擇雙尾檢驗[20];主成分分析選擇KMO和Bartlett法進行球形度檢驗,利用最大方差法獲得旋轉因子,旋轉方法為Kaiser標準化的正交旋轉法,旋轉在迭代5次后收斂[21-22];聚類分析選擇按照組間連接系統(tǒng)聚類,選擇Pearson相關性作為度量標準.
關聯(lián)性顯著的元素可能具有同源性,能夠提高預測精度[23-24],但是即使部分污染物關聯(lián)性較弱,二者之間可能存在非線性相關性,仍可能對預測精度有積極影響. 因此,該文中關聯(lián)性較弱的污染物數(shù)據(jù)仍作為神經(jīng)網(wǎng)絡預測模型的輸入因子進行訓練,并設置將關聯(lián)性較弱因子去除后的輸入因子預測模型作為對照樣本,探究關聯(lián)性弱的因子對BP神經(jīng)網(wǎng)絡模型預測精度的影響.
BP神經(jīng)網(wǎng)絡是一種典型的多層前饋網(wǎng)絡[25],其函數(shù)逼近、回歸計算等能力已被廣泛應用于環(huán)境科學研究領域,在預測土壤污染物的空間分布方面均取得了較好的效果[26-28].
該研究所選場地范圍較大,分析樣品數(shù)量多,受成本限制,個別采樣點的某些污染物指標未檢測,導致數(shù)據(jù)缺失. 為保證檢測數(shù)據(jù)的完整性,通過構建BP神經(jīng)網(wǎng)絡對缺失數(shù)據(jù)進行預測:將樣本中不含缺失數(shù)據(jù)的因子作為模型的輸入條件,含有缺失數(shù)據(jù)的因子作為輸出條件,利用樣本中的已知數(shù)據(jù)訓練神經(jīng)網(wǎng)絡,在訓練達到要求后,將缺失樣本的已知數(shù)據(jù)輸入模型,輸出值即為缺失數(shù)據(jù)的預測值[29]. 研究區(qū)共獲得1 691組土壤樣本數(shù)據(jù),隨機抽取30組作為驗證樣本,在其余樣本中隨機選取50、200、800和1 661組(剩余全部樣本)作為訓練樣本,其中,w(Cu)、w(Pb)和w(DBA)、w(BkF)、w(BbF)、w(BaA)、w(Nap)、w(Chr)作為缺失數(shù)據(jù),通過構建的神經(jīng)網(wǎng)絡對缺失數(shù)據(jù)進行預測,并與驗證樣本的實測數(shù)據(jù)進行對比,評價模型的預測精度.
重金屬含量預測選取采樣點空間坐標參數(shù)、土壤pH、含水率(ω)、w(As)、w(Cr)、w(Zn)、w(Ni)、w(Cd)作為輸入?yún)?shù),w(Pb)和w(Cd)作為輸出參數(shù);PAHs含量預測選取采樣點空間坐標參數(shù)、土壤pH、含水率(ω)和w(BaP)作為輸入?yún)?shù),w(DBA)、w(BkF)、w(BbF)、w(BaA)、w(Nap)、w(Chr)作為輸出參數(shù).
構建的BP神經(jīng)網(wǎng)絡拓撲結構見圖2,神經(jīng)網(wǎng)絡每個節(jié)點代表一種特定的輸出函數(shù),稱為激活函數(shù)[30],每兩個節(jié)點的連接代表權重,通過不斷學習來調整權重值[31],該文采用的學習速率為0.05. 構建的雙層隱含層單元數(shù)分別為20和15個. 選取正切函數(shù)tansig作為隱含層的傳遞函數(shù),線性函數(shù)purelin作為輸出層的傳遞函數(shù),共軛梯度函數(shù)trainscg作為訓練函數(shù)對樣本進行訓練[32]. 重金屬和PAHs指標最大迭代次數(shù)分別設置為20 000和 10 000 次,通過反復迭代,最終確定權值和閾值,建立預測模型.
圖2 BP神經(jīng)網(wǎng)絡拓撲結構Fig.2 Topological of BP neural network
BP神經(jīng)網(wǎng)絡根據(jù)梯度下降法調節(jié)連接權值,使訓練誤差達到最小,為使輸入數(shù)據(jù)能夠適應傳遞函數(shù),提高計算過程的收斂速度,需要對訓練數(shù)據(jù)進行歸一化處理[33],處理后的變量取值范圍在[-1,1]之間,計算方法如式(1)所示:
(1)
式中,xi為輸入變量ti歸一化后的值,tmax和tmin分別為輸入變量X的最大值和最小值.
所有試驗數(shù)據(jù)使用Microsoft Excel 2013進行處理,并采用SPSS 22軟件進行多元統(tǒng)計分析;利用Matlab 2014軟件編寫程序構建BP神經(jīng)網(wǎng)絡的預測模型,采用OriginPro 9.0軟件進行模型精度分析.
研究區(qū)土壤污染物含量統(tǒng)計結果見表1[34]. 結果顯示,除w(Chr)外,土壤中重金屬和PAHs均存在不同程度的超標現(xiàn)象. 土壤重金屬污染物中,除w(Cr)外,其他6種重金屬平均含量均高于當?shù)赝寥辣尘爸?,與GB 36600—2018《土壤環(huán)境質量 建設用地土壤污染風險管控標準》[35]中的風險篩選值對比,7種重金屬含量超標率表現(xiàn)為w(Ni)>w(Cu)>w(As)>w(Pb)>w(Zn)=w(Cd)>w(Cr),最大超標倍數(shù)分別為43.93、8.45、5.80、8.05、5.03、17.70和2.98倍;土壤PAHs污染物中,w(BaP)的超標情況較為嚴重,局部超標倍數(shù)較高,與GB 36600—2018中土壤PAHs風險篩選值對比,除w(Chr)未超標外,其他6種PAHs含量超標率表現(xiàn)為w(BaP)>w(Chr)>w(DBA)>w(BbF)=w(BaA)>w(Nap)>w(BkF),最大超標倍數(shù)分別為120.00、28.00、30.55、30.36、3.02和1.22倍.
變異系數(shù)(CV)可以反映土壤污染物的變異程度,該值越大,表明該污染物受人為活動的影響越強[36]. 研究區(qū)土壤重金屬和PAHs的變異系數(shù)見表1[34],其中重金屬Zn、Cu、Pb、Ni和Cd的變異系數(shù)范圍為2.75~8.12,屬于強變異(CV≥1.0);As和Cr的變異系數(shù)分別為0.43和0.82,屬于中等變異(0.1≤CV<1.0);PAHs變異系數(shù)為3.49~6.84,屬于強變異[37]. 由此可見,場地重金屬和PAHs污染物受到較強的人為污染源影響.
表1 研究區(qū)土壤重金屬和PAHs的含量統(tǒng)計值[34]
2.2.1重金屬污染物
相關性分析通常作為污染物同源鑒別的依據(jù)[38-39]. 從表2可以得出,Zn與Cu、Pb、Cd均具有較強相關性(R為0.579~0.673,P<0.01),Cu與Ni相關性較好(R=0.519),Cr與其他6種重金屬相關性均較差,其他重金屬間彼此相關性較弱,表明Cr與其他重金屬的來源均不同. 由于研究場地不同類型的生產車間較多,因此不同區(qū)域土壤重金屬污染來源差異性較大.
表2 研究區(qū)土壤7種重金屬含量的相關性
KMO檢驗值(0.607)和Bartlett球形度(ξ=4 045,P<0.05)分析結果表明,該研究中重金屬污染物含量數(shù)據(jù)適用主成分分析. 7種重金屬經(jīng)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,簇距離越低,表示因子間關聯(lián)越顯著[40-41]. 由圖3(b)可見,第一簇中Zn與Pb、As、Cd的簇距離較小,表明Zn與這3種重金屬具有相似污染源;Cr與其他簇的距離均較大,表明Cr與其他6種重金屬的來源差異性較大.
圖3 研究區(qū)土壤中7種重金屬主成分和聚類分析Fig.3 Principal component diagram and dendrogram of seven heavy metals contents in soils of the study area
上述3種多元統(tǒng)計分析方法所得結果相近,綜合分析結果并結合場地生產歷史可以得出:不同生產車間的特征污染物不同,Zn、Pb、As和Cd是熔鑄車間的主要產品或原料,關聯(lián)性較好;Cu和Ni主要是管棒車間的主要產物,關聯(lián)性較好,這些重金屬污染主要源自不同生產線的原材料和產品;Cr僅在盤管車間涉及污染,且w(Cr)(平均值為54.9 mg/kg)與當?shù)赝寥辣尘爸?57.7 mg/kg)接近,相關研究表明,Cr污染主要源于成土母質[42-43],研究區(qū)Cr判斷為受少量工業(yè)影響的自然源,因此與其他6種重金屬關聯(lián)性較弱.
2.2.2PAHs污染物
研究區(qū)土壤中7種PAHs相關性分析結果見表3,除Nap外,BaP、DBA、BkF、BbF、BaA和Chr彼此間均呈顯著相關(P<0.01),表明這6種PAHs具有相同或相似的來源;Nap與其他PAHs相關性較弱.
表3 研究區(qū)土壤中7種PAHs的相關性
KMO檢驗值(0.854)和Bartlett球形度(ξ=23 761,P<0.05)分析表明,PAHs含量數(shù)據(jù)適用主成分分析. 通過分析獲得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 研究區(qū)土壤中7種PAHs主成分和聚類分析Fig.4 Principal component diagram and dendrogram of seven PAHs in soils of the study area
結合廠區(qū)生產歷史并綜合3種分析方法結果表明,除Nap外,BaP、DBA、BkF、BbF、BaA和Chr彼此關聯(lián)性較好,其中BaP、BaA和Chr是燃煤型排放的標志污染物,BkF、BbF和DBA主要來自汽油或柴油的燃燒[44],因此推斷這6種PAHs主要來自燃煤和交通混合源[45];而Nap的產生可能與生物質不完全燃燒有關[46],因此與其他6種PAHs關聯(lián)性較弱.
2.3.1不同訓練樣本數(shù)對神經(jīng)網(wǎng)絡精度影響
按照1.4節(jié)建立的BP神經(jīng)網(wǎng)絡模型進行訓練,得到不同訓練樣本數(shù)量對訓練精度的影響如圖5所示. 由圖5可以看出,隨著訓練樣本數(shù)量的增加,神經(jīng)網(wǎng)絡訓練誤差(mean squared error,MSE)逐漸減小,并且訓練結束時MSE更接近設置的目標誤差,表明隨著樣本數(shù)量的增加,訓練精度逐步提高;而當訓練樣本數(shù)據(jù)減至50個時,重金屬和PAHs相應BP模型在未達到迭代次數(shù)(20 000 和 10 000 次)時提前終止訓練,表明隨著迭代次數(shù)的增加,訓練誤差不再降低,訓練精度無法提高,停止訓練. 因此,適當增加訓練樣本數(shù)量能夠有效提高模型預測精度.
圖5 不同訓練樣本數(shù)量對BP神經(jīng)網(wǎng)絡精度的影響Fig.5 The effect of different training sample size on the accuracy of BP neural network
2.3.2神經(jīng)網(wǎng)絡訓練精度和誤差分析
選擇全部1 661個訓練樣本建立BP神經(jīng)網(wǎng)絡模型,結果如圖6所示. 重金屬及PAHs樣本在訓練次數(shù)分別為 20 000 和 10 000 次時達到最佳效果,訓練樣本MSE分別為1.9×10-4和4.3×10-5,表明訓練收斂效果較好;同時,訓練樣本輸出值和目標值之間的相關系數(shù)(R)分別達到 0.990 11 和 0.995 19,R接近1,表明重金屬和PAHs樣本神經(jīng)網(wǎng)絡訓練效果較好.
圖6 BP神經(jīng)網(wǎng)絡訓練誤差曲線和線性回歸結果Fig.6 The error curve diagram and linear regression diagram of BP neural network training
利用訓練完成后的BP神經(jīng)網(wǎng)絡模型對30組驗證樣本的缺失數(shù)據(jù)〔w(Cu)、w(Pb)和w(DBA)、w(BkF)、w(BbF)、w(BaA)、w(Nap)、w(Chr)〕進行預測,結果如圖7所示. 缺失數(shù)據(jù)預測結果與實測結果的R范圍為0.901~0.996,預測曲線與實測曲線重合度較好. 將BP神經(jīng)網(wǎng)絡模型預測結果與實測結果進行對比分析,結果(見表4)顯示,驗證樣本的決定系數(shù)(R2)范圍為0.812~0.993,表明該模型預測擬合度效果較高. 模擬效率系數(shù)(Nash-Sutcliffe efficiency coefficient,NSE)通常是衡量模擬結果可靠性的重要指標[29],取值范圍為(-∞,1],當NSE為正值時,表明模擬結果可信,且該值越接近1,表示模型匹配程度越好. 由表4可見:該研究中NSE范圍為0.779~0.959,除w(DBA)外,其他驗證樣本的NSE均大于0.8,模擬效果較好,這與其他相關研究結果[47]相似;同時,驗證樣本的均方根誤差(RMSE)和平均絕對誤差(MAE)均較小. 綜上,BP神經(jīng)網(wǎng)絡模型可靠性高,能夠較好地預測研究區(qū)土壤重金屬和PAHs含量.
圖7 BP神經(jīng)網(wǎng)絡預測值與實測值結果散點圖Fig.7 The scatter diagram of predicted values and measured values of BP neural network
表4 BP神經(jīng)網(wǎng)絡預測精度和誤差分析
為了探究污染物關聯(lián)性對BP神經(jīng)網(wǎng)絡模型預測精度的影響,根據(jù)污染物的關聯(lián)性分析結果,該研究去除了關聯(lián)性較弱的因子〔w(Cd)、w(Cr)〕,選取重金屬中關聯(lián)性較強的因子〔(w(As)、w(Cr)、w(Zn)、w(Ni)〕和采樣點空間參數(shù)以及理化性質作為輸入變量,w(Cu)、w(Pb)作為輸出變量,在其他參數(shù)不變的條件下構建對照模型,該對照模型的預測結果與選取全部因子作為輸入變量構建的模型預測結果對比見表5. 由表5可以看出,與去除關聯(lián)性較弱的因子相比,選取全部因子作為輸入變量建立的預測模型整體擬合度更好,其精度分析指標(R2)相比于對照模型提高了0.129,NSE提高了0.134,誤差分析指標RMSE和MAE分別降低了0.036和0.031. 由此可見,相比于僅使用關聯(lián)性強的因子作為輸入?yún)?shù),考慮不同影響因子構建的BP神經(jīng)網(wǎng)絡模型預測更加準確可靠. 推測原因可能是,關聯(lián)性較弱因子之間仍有非線性相關性,但這種關聯(lián)難以通過具體數(shù)學方程計算,而神經(jīng)網(wǎng)絡不需要確定具體的映射方程,因此將這些因子作為輸入變量構建神經(jīng)網(wǎng)絡,能夠和輸出端建立較好的映射關系,從而進一步提高預測模型的精度.
表5 不同輸入變量的預測模型精度對比
a) 通過分析某金屬加工廠生產區(qū)內土壤污染特性,結果表明,除PAHs的Chr外,重金屬和PAHs均存在不同程度的超標現(xiàn)象,其中重金屬Ni、As、Cu及PAHs中BaP超標較為嚴重;重金屬和PAHs變異系數(shù)較高,該區(qū)域土壤受到較強的人為污染源影響.
b) 多元統(tǒng)計分析結果顯示,污染土壤中重金屬Zn與Pb、As、Cd關聯(lián)性均較好,Cu與Ni關聯(lián)性較好,而Cr與其他6種重金屬關聯(lián)性較弱,Zn、Pb、As、Cd、Cu和Ni污染主要源于不同生產線的原材料和產品,Cr為受少量工業(yè)影響的自然源;PAHs中除Nap外,BaP、DBA、BkF、BbF、BaA和Chr彼此關聯(lián)性較好,6種關聯(lián)性好的PAHs主要來自燃煤和交通混合源,Nap則源于生物質不完全燃燒.
c) 構建的BP神經(jīng)網(wǎng)絡模型訓練效果較好;驗證樣本的預測值與實測值之間的誤差分析顯示,各污染物含量的決定系數(shù)(R2)范圍為0.812~0.993,模擬效率系數(shù)(NSE)范圍為0.779~0.959,均方根誤差(RMSE)和平均絕對誤差(MAE)均較?。粯嫿ǖ腂P神經(jīng)網(wǎng)絡模型準確可靠,能夠較好地預測研究區(qū)內土壤污染物含量,且關聯(lián)性較弱因子的輸入能進一步提高預測模型的精度.