龍口市污灌區(qū)重金屬污染影響因子分析*

劉建華 顧振飛 曹見飛，2 徐明雪 王翠秀# 吳泉源

(1.山東師范大學(xué)地理與環(huán)境學(xué)院，山東 濟(jì)南 250358；2.偉志股份公司，福建 泉州 362200)

面對(duì)全球水資源匱乏的現(xiàn)狀，污水灌溉已經(jīng)成為許多地區(qū)農(nóng)田灌溉的主要方式，然而污水中含有的重金屬元素將不斷在土壤中富集，并會(huì)通過農(nóng)作物進(jìn)入食物鏈影響人類健康，由此產(chǎn)生的風(fēng)險(xiǎn)引起學(xué)者們的廣泛關(guān)注[1-2]。對(duì)污灌區(qū)土壤重金屬進(jìn)行影響因子分析有助于重金屬污染的有效治理，對(duì)預(yù)防重金屬污染產(chǎn)生的潛在危害也具有重要的意義。污灌區(qū)土壤重金屬的含量同時(shí)受到多個(gè)因子的影響，傳統(tǒng)的影響因子研究方法主要有Pearson相關(guān)系數(shù)法、主成分分析法、地理加權(quán)回歸等方法[3-5]。這類回歸統(tǒng)計(jì)方法雖然可以測(cè)量和估計(jì)各影響因子的影響程度，但這些分析方法都是在假設(shè)預(yù)測(cè)變量獨(dú)立的情況下進(jìn)行的，忽略了驅(qū)動(dòng)因素和重金屬濃度之間的空間關(guān)系，并且受影響因子間多重共線性的阻礙，從而產(chǎn)生較大的方差，降低了模型的推理能力[6]。2010年，WANG等[7]研發(fā)了地理探測(cè)器模型并用于地理環(huán)境因子分析，目前已有多名學(xué)者將地理探測(cè)器模型引入土壤重金屬污染的影響因子解析中，并且取得理想的效果[8-11]。為此，本研究采用地理探測(cè)器模型綜合分析污灌區(qū)重金屬污染的影響因子，解決以往研究中影響因子解釋度低的問題。

龍口市是山東省重要的老工業(yè)基地，同時(shí)也是膠東農(nóng)產(chǎn)品生產(chǎn)基地，該地污水灌溉歷史悠久，污水中的有害物質(zhì)在土壤中不斷積累，使得農(nóng)田生態(tài)安全面臨嚴(yán)峻的挑戰(zhàn)，農(nóng)田污染現(xiàn)狀調(diào)查及農(nóng)田生態(tài)影響分析等都是亟待解決的問題。本研究將重金屬富集因子特征分析與地理探測(cè)器模型分析相結(jié)合，評(píng)估龍口市污灌區(qū)重金屬的污染水平與空間分布，分析各影響因子對(duì)重金屬污染的貢獻(xiàn)，探索不同影響因子之間的協(xié)同或拮抗作用，在此基礎(chǔ)上，比較了地理探測(cè)器模型與Pearson相關(guān)系數(shù)法的優(yōu)缺點(diǎn)，以期為污灌區(qū)重金屬空間分布特征與影響因子分析提供一種高速有效的方法，為農(nóng)田規(guī)劃與環(huán)境治理提供借鑒。

1 材料與方法

1.1 研究區(qū)概況

研究區(qū)位于山東省龍口市北部平原地區(qū)(37°34′35″N～37°44′12″N， 120°20′56″E～120°40′47″E)，區(qū)域總面積409 km2。研究區(qū)內(nèi)成土母質(zhì)以沖積和海積物為主，土壤類型主要包括褐土、棕壤和風(fēng)砂，pH為6.5～8.0，成土條件基本一致。研究區(qū)距海較近，屬溫帶季風(fēng)性氣候，年平均氣溫12 ℃左右，年均降水量600 mm左右，年降水分配不均，主要集中在5—9月，具有春秋季旱、夏季澇的特點(diǎn)[12]。研究區(qū)內(nèi)大部分區(qū)域?yàn)檗r(nóng)田，主要種植小麥、玉米等作物，核心建成區(qū)主要分布著蘆頭工業(yè)區(qū)、住宅區(qū)、商業(yè)區(qū)等。研究區(qū)內(nèi)煤炭、金、高嶺土、砂礦等礦產(chǎn)資源豐富，煤炭、采礦、鑄造和電鍍等工業(yè)發(fā)達(dá)，其中工業(yè)殘?jiān)亩逊?、廢液排放都會(huì)對(duì)周圍的環(huán)境產(chǎn)生一定的影響。由于當(dāng)?shù)亟邓y以滿足農(nóng)業(yè)需求，經(jīng)污水處理廠處理后的生活廢水、工業(yè)污水經(jīng)由黃水河及永汶河后進(jìn)入各支渠用于農(nóng)田灌溉，加劇農(nóng)田土壤中各類重金屬元素的富集，區(qū)域內(nèi)主要污染重金屬為Cr、Ni、Pb、Zn、Cu、Co、As、Cd。

1.2 樣品采集和理化分析

于2020年5、11月采集土壤樣品，采樣期間降水稀少，且農(nóng)作物處于成長(zhǎng)的關(guān)鍵期，一般會(huì)采用河水灌溉。根據(jù)研究區(qū)各種土地利用類型占比及主要污染源(大型工礦業(yè))分布情況布設(shè)采樣點(diǎn)，樣品采集深度為0～20 cm，以30 m×30 m為一個(gè)樣方，清理土壤表面后用小木鏟進(jìn)行采集，四分法混合均勻后剩余約1 kg裝入聚乙烯自封袋中低溫封裝帶回。本研究共設(shè)70個(gè)采樣點(diǎn)，皆采用全球定位系統(tǒng)(GPS)定位，采樣點(diǎn)分布見圖1。首先挑除過濾土壤樣品中的粗大雜質(zhì)，然后使用瑪瑙研缽磨細(xì)直至通過200目孔徑篩。對(duì)預(yù)處理后的土壤樣品進(jìn)行重金屬檢測(cè)分析，采用石墨爐原子吸收光譜法(GF-AAS)測(cè)定Cd含量；用電感耦合等離子體原子發(fā)射光譜法(ICP-OES)測(cè)定Cu、Pb、Zn、Co、Ni含量；采用氫化物發(fā)生-原子熒光光譜法(HG-AFS)測(cè)定As含量；采用火焰原子吸收分光光度法測(cè)定Cr、Mn含量。每個(gè)樣品重復(fù)測(cè)量3次取平均值，用去離子水調(diào)節(jié)土壤樣品土水比為1 g∶5 mL后測(cè)定pH。最后，根據(jù)國(guó)家標(biāo)準(zhǔn)土壤樣品(GBW7401)進(jìn)行準(zhǔn)確度和精密度的控制，測(cè)量結(jié)果均符合質(zhì)量控制要求。

圖1 研究區(qū)采樣點(diǎn)分布Fig.1 Distribution of sampling points in the study area

1.3 數(shù)據(jù)來源與統(tǒng)計(jì)分析

研究區(qū)位于龍口市污灌區(qū)，灌溉系統(tǒng)和農(nóng)業(yè)投入品對(duì)當(dāng)?shù)赝寥乐亟饘傥廴揪嬖谝欢ㄓ绊?。但研究區(qū)范圍較小，農(nóng)產(chǎn)品相對(duì)單一，各采樣點(diǎn)灌溉習(xí)慣和農(nóng)業(yè)投入情況差異較小，因此不將灌溉系統(tǒng)和農(nóng)業(yè)投入品作為影響因子進(jìn)行研究。本研究主要考察7個(gè)影響因子對(duì)重金屬污染的影響，其中包括3個(gè)人為因素(工業(yè)密度(DI)、人口密度(DP)和道路密度(DR))和4個(gè)自然因素(植被覆蓋度(NDVI)、土壤母質(zhì)(PMT)、土壤有機(jī)質(zhì)(SOM)和高程(DEM))。其中，NDVI、PMT、DEM數(shù)據(jù)來源于2014年山東省土壤類型分布圖；DP數(shù)據(jù)來源于龍口市統(tǒng)計(jì)年鑒；使用ArcGIS 10.0軟件采樣點(diǎn)為中心建立2 000 m的緩沖區(qū)，計(jì)算每個(gè)緩沖區(qū)內(nèi)的DI和DR數(shù)據(jù)；SOM為點(diǎn)數(shù)據(jù)，與70個(gè)土壤采樣點(diǎn)一一對(duì)應(yīng)，實(shí)測(cè)獲得。采用地理探測(cè)器模型進(jìn)行分析時(shí)需要先將數(shù)據(jù)離散化，因此采用ArcGIS 10.0軟件中的自然斷裂分類方法將連續(xù)數(shù)據(jù)轉(zhuǎn)換為離散數(shù)據(jù)，以滿足地理探測(cè)器模型的數(shù)據(jù)格式要求。為了評(píng)價(jià)地理探測(cè)器模型的效果，將分析結(jié)果與Pearson相關(guān)系數(shù)法的分析結(jié)果進(jìn)行比較。

1.4 研究方法

1.4.1 富集系數(shù)

富集系數(shù)可以用來粗略評(píng)價(jià)環(huán)境中重金屬的污染程度(富集程度)并分析重金屬來源，富集系數(shù)計(jì)算見式(1):

(1)

式中:EFi為重金屬i的富集系數(shù)；Ci為土壤樣品中重金屬i的質(zhì)量濃度，mg/kg；Cn為土壤樣品中參考元素的質(zhì)量濃度，mg/kg；Bi為重金屬i的背景值，mg/kg；Bn為參考元素的背景值，mg/kg。

對(duì)研究區(qū)各重金屬元素的含量特征進(jìn)行分析發(fā)現(xiàn)，Mn的變異性較小，受人為活動(dòng)影響較小[13]，且含量水平基數(shù)與其他元素相比較大，因此本研究中選擇Mn作為參考元素。在重金屬來源分析時(shí)，認(rèn)為富集系數(shù)≤1.5時(shí)重金屬主要受自然因素影響，當(dāng)富集系數(shù)>1.5時(shí)重金屬主要受人為因素影響[14]。

1.4.2 地理探測(cè)器模型

本研究采用地理探測(cè)器模型分析影響因子對(duì)研究區(qū)重金屬污染的影響。地理探測(cè)器模型是一種分析地理事物空間分異及其驅(qū)動(dòng)力的空間統(tǒng)計(jì)方法，可用于分析因子的交互作用對(duì)因變量的影響。地理探測(cè)器模型由風(fēng)險(xiǎn)探測(cè)器、因子探測(cè)器、生態(tài)探測(cè)器和交互探測(cè)器4大模塊組成，通過計(jì)算單個(gè)或多個(gè)因子與因變量間的解釋力(q)分析因子對(duì)因變量的影響程度。q的取值范圍為0～1，q=0時(shí)表示因子與因變量間無關(guān)聯(lián)，q=1時(shí)表示因變量可全部由該因子解釋，q越大表明因子對(duì)因變量的影響越大。在因子交互探測(cè)中，交互探測(cè)模塊通過比較2個(gè)單因子a1、a2交互后的解釋力(即q(a1∩a2))與a1、a2單獨(dú)作用時(shí)的解釋力(即q(a1)、q(a2))大小，判別兩者間是否存在交互作用及其作用強(qiáng)弱，交互探測(cè)表達(dá)式見表1[15]。

表1 交互探測(cè)表達(dá)式Table 1 Expression of interactive detection

2 結(jié)果與討論

2.1 污灌區(qū)重金屬污染的描述性統(tǒng)計(jì)與評(píng)價(jià)

研究區(qū)土壤樣品pH為7.00～7.85，平均值為7.43，屬于中性偏堿性土壤，土壤樣品中重金屬質(zhì)量濃度與富集系數(shù)統(tǒng)計(jì)見表2。Cu、Cr、Ni、Pb、Zn、As、Cd、Co的平均值均高于龍口市土壤背景值，最大值分別為背景值的9.87、1.92、3.17、7.70、3.35、3.25、8.06、1.79倍，說明研究區(qū)土壤中重金屬受人為干擾較大，集聚趨勢(shì)比較明顯。本研究選取Mn作為參考元素，計(jì)算得到研究區(qū)Cu、Cr、Ni、Pb、Zn、As、Cd、Co的富集系數(shù)平均值分別1.96、1.32、1.49、1.75、1.69、1.56、3.76、1.28。其中，Cu、Pb、Zn、As、Cd的富集系數(shù)大于1.5，可以被認(rèn)為是受到人為因素影響的主要重金屬元素，因此在研究區(qū)中富集現(xiàn)象明顯；Co、Cr、Ni富集系數(shù)低于1.5，且在研究區(qū)范圍內(nèi)標(biāo)準(zhǔn)偏差均小于等于0.2，說明這3種重金屬的富集主要受自然因素的影響較大，受到人為因素影響較小。

表2 研究區(qū)重金屬質(zhì)量濃度和富集因子統(tǒng)計(jì)Table 2 Statistical of mass concentration of heavy metals and enrichment factor in study area

2.2 重金屬富集系數(shù)的空間分布

圖2為研究區(qū)70個(gè)采樣點(diǎn)8種重金屬富集系數(shù)的空間分布情況。8種重金屬的富集系數(shù)分布趨勢(shì)總體相似，基本表現(xiàn)為西北部富集系數(shù)相對(duì)偏高，這是因?yàn)樵摰貐^(qū)位于東海工業(yè)園區(qū)附近(發(fā)電廠東部)；此外，富集系數(shù)的高值點(diǎn)位于研究區(qū)中東部，靠近徐福鎮(zhèn)的煤炭開采地區(qū)域。8種重金屬中，Cd的富集系數(shù)普遍偏高，表明該元素污染較為嚴(yán)重，可能是農(nóng)資產(chǎn)品、地膜以及除草劑在農(nóng)田中普遍使用造成的[16-17]；Co、Cr的富集系數(shù)在全區(qū)范圍內(nèi)均偏小，富集程度較輕?？傮w看來，研究區(qū)重金屬富集系數(shù)的分布情況與工業(yè)區(qū)的分布具有相似性，主要受到工業(yè)園區(qū)、煤礦區(qū)、農(nóng)作種植等的影響。

2.3 地理探測(cè)器模型分析

2.3.1 因子分析

利用地理探測(cè)器模型中的因子探測(cè)器模塊評(píng)估7個(gè)影響因子對(duì)8種重金屬污染的解釋力，結(jié)果見表3?？傮w看來，7種影響因子中DI對(duì)8種重金屬的平均解釋力最大,其次為SOM和PMT，說明工業(yè)是重金屬污染的主要來源。實(shí)地調(diào)研發(fā)現(xiàn)，龍口市許多工廠、煤礦區(qū)過去直接向河流系統(tǒng)排放工業(yè)廢水，以往流域內(nèi)的發(fā)電廠、電鍍、冶煉等企業(yè)也存在將未經(jīng)處理或處理不徹底的污水排放至河流及農(nóng)田的情況，因此研究區(qū)內(nèi)工業(yè)分布是造成重金屬污染的主要原因[18]。8種重金屬元素中，DI對(duì)Cd、Pb的解釋力較大，說明工業(yè)是影響土壤樣品Cd、Pb污染的主要影響因子。研究區(qū)內(nèi)煤礦區(qū)分布較多，長(zhǎng)期以來的煤炭開采、煤矸石的堆放使污染物隨雨水進(jìn)入河流，煤炭廢渣以及汽車尾氣排放的Cd、Pb經(jīng)過大氣的干濕沉降遷入土壤，這些都是造成附近土壤中Pb、Cd富集的主要原因[19]；此外，農(nóng)業(yè)活動(dòng)中為提高產(chǎn)量而使用的尿素、磷肥也會(huì)使得土壤中Pb、Cd含量增加[20-21]。SOM對(duì)Cu和Zn的解釋力較大，主要是因?yàn)镾OM對(duì)重金屬有較強(qiáng)的絡(luò)合作用。有研究表明，豬飼料中添加了大量的Cu、Zn來防治疾病、加快生長(zhǎng)，但是90%(質(zhì)量分?jǐn)?shù))以上會(huì)以豬糞的形式排出，且豬糞中有機(jī)質(zhì)豐富，絡(luò)合作用較強(qiáng)[22-23]，因此SOM是土壤中Cu、Zn富集的主要影響因子。Cr受到PMT、DP、SOM及DI共同作用，說明Cr污染的影響因子較為復(fù)雜；Co、As和Ni等主要受PMT的影響，人類活動(dòng)對(duì)這類重金屬在土壤中的富集影響較小，這與前人的研究結(jié)果總體一致[24]。本研究中DR對(duì)各重金屬富集的影響處于較低水平，這可能與大多數(shù)工業(yè)區(qū)位于郊區(qū)，DR通常較低，這在某種程度上可能掩蓋了DR對(duì)重金屬污染的影響。

圖2 研究區(qū)8種重金屬富集系數(shù)的空間分布Fig.2 Spatial distribution of 8 heavy metals enrichment factor in study area

2.3.2 交互作用分析

從圖2(g)可以看出，Pb分布的空間異質(zhì)性較大，在影響因子分析中更具有代表性，因此以Pb為例考察影響因子間對(duì)重金屬污染的交互作用。表4為影響因子間交互作用對(duì)Pb污染的解釋力，通過比較單一因子解釋力和交互作用解釋力的大小關(guān)系，判斷兩個(gè)影響因子的交互作用效果。經(jīng)比較，NDVI、SOM、DI、DR、DP、DEM、PMT間的交互作用均為非線性增強(qiáng)，也即影響因子間單因子解釋力合計(jì)小于影響因子交互作用的解釋力，Pb分布受到各影響因子間交互作用的影響。其中，DI和DP交互作用對(duì)Pb的解釋力達(dá)0.921，是影響Pb分布的主要因子，其他影響因子也會(huì)加強(qiáng)工業(yè)以及人口因素對(duì)Pb分布的影響。

影響因子間交互作用非線性增強(qiáng)，顯示了人為因素和自然因素在控制污灌區(qū)重金屬污染上存在相互作用而非相互獨(dú)立，影響因子間可能通過復(fù)雜的聯(lián)動(dòng)機(jī)制增強(qiáng)對(duì)污灌區(qū)重金屬濃度的影響。對(duì)其他重金屬的影響因子分析時(shí)也發(fā)現(xiàn)了相同的規(guī)律，可見研究區(qū)重金屬污染是多個(gè)影響因子共同作用的結(jié)果。

2.4 地理探測(cè)器模型與Pearson相關(guān)系數(shù)的比較

Pearson相關(guān)系數(shù)法是影響因子相關(guān)性分析的常用方法，因此采用Pearson相關(guān)系數(shù)法來驗(yàn)證地理探測(cè)器模型分析結(jié)果的可靠性，并以Pb影響因子分析為例進(jìn)行對(duì)比。由圖3可見，兩種方法在確定最主要影響因子時(shí)結(jié)果一致，均認(rèn)為DI是Pb污染的主要影響因子，其次為NDVI。然而，NDVI、DR、DP、PMT的Pearson相關(guān)系數(shù)為負(fù)，但地理探測(cè)器模型得到的解釋力為正。這是因?yàn)榈乩硖綔y(cè)器模型對(duì)影響因子相互關(guān)系的探測(cè)既包括線性關(guān)系也存在非線性關(guān)系，而Pearson相關(guān)系數(shù)探討的是影響因子與重金屬間的線性關(guān)系，分為正相關(guān)和負(fù)相關(guān)，如果影響因子與重金屬之間出現(xiàn)非線性關(guān)系時(shí)，則會(huì)出現(xiàn)錯(cuò)誤的判斷。總體而言，Pearson相關(guān)系數(shù)法驗(yàn)證了地理探測(cè)器模型對(duì)污灌區(qū)中Pb污染主要影響因子的識(shí)別，但在較弱影響因子評(píng)估中存在一定差異。在對(duì)Zn、Pb、Cu、Cr、Ni、As的影響因子分析中，地理探測(cè)器模型和Pearson相關(guān)系數(shù)法的評(píng)估結(jié)果幾乎相同，但對(duì)Co、Cd的影響因子分析中有所偏差。與Pearson相關(guān)系數(shù)法相比，地理探測(cè)器模型在檢測(cè)重金屬污染和潛在影響因子間的非線性關(guān)系中具有一定優(yōu)勢(shì)。

表3 各影響因子對(duì)8種重金屬污染的解釋力Table 3 The explanatory power of different influence factors for 8 heavy metals pollution

表4 影響因子間交互作用對(duì)Pb污染的解釋力Table 4 The explanatory power of interactions between influence factors for Pb pollution

圖3 地理探測(cè)器模型與Pearson相關(guān)系數(shù)法的比較Fig.3 Comparison of geographical detector model and Pearson correlation index method

3 結(jié) 論

(1) 研究區(qū)內(nèi)土壤重金屬呈現(xiàn)不同程度的富集，從富集系數(shù)平均值來看，Cu、Pb、Zn、As、Cd的富集程度較大，主要受到人為因素的影響；Co、Cr、Ni富集程度較低，主要受自然因素的影響較大，受人為因素影響較小。從富集系數(shù)空間分布來看，8種重金屬富集系數(shù)分布情況相似，主要受到工業(yè)園區(qū)、煤礦區(qū)、農(nóng)作種植等的影響。

(2) 地理探測(cè)器因子探測(cè)結(jié)果表明，7種影響因子中DI對(duì)8種重金屬的平均解釋力最大,其次為SOM和PMT。不同重金屬的主要影響因子也存在不同，Cd、Pb主要?dú)w因于DI，Zn、Cu主要?dú)w因于SOM，As、Co、Ni主要?dú)w因于PMT，Cr分布受到PMT、DP、SOM、DI的共同作用影響。

(3) 基于地理探測(cè)器模型交互作用分析，不同影響因子之間的交互作用強(qiáng)弱不同但均為非線性增強(qiáng)，污灌區(qū)重金屬污染是多種影響因子共同作用的結(jié)果。

(4) 與Pearson相關(guān)系數(shù)法比，地理探測(cè)器模型可以有效分析因子間的交互作用對(duì)重金屬污染的影響，在檢測(cè)潛在影響因子間的非線性關(guān)系中具有一定優(yōu)勢(shì)。