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

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

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

龍口市是山東省重要的老工業(yè)基地，同時也是膠東農(nóng)產(chǎn)品生產(chǎn)基地，該地污水灌溉歷史悠久，污水中的有害物質(zhì)在土壤中不斷積累，使得農(nóng)田生態(tài)安全面臨嚴峻的挑戰(zhàn)，農(nóng)田污染現(xiàn)狀調(diào)查及農(nóng)田生態(tài)影響分析等都是亟待解決的問題。本研究將重金屬富集因子特征分析與地理探測器模型分析相結(jié)合，評估龍口市污灌區(qū)重金屬的污染水平與空間分布，分析各影響因子對重金屬污染的貢獻，探索不同影響因子之間的協(xié)同或拮抗作用，在此基礎(chǔ)上，比較了地理探測器模型與Pearson相關(guān)系數(shù)法的優(yōu)缺點，以期為污灌區(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ì)以沖積和海積物為主，土壤類型主要包括褐土、棕壤和風砂，pH為6.5～8.0，成土條件基本一致。研究區(qū)距海較近，屬溫帶季風性氣候，年平均氣溫12 ℃左右，年均降水量600 mm左右，年降水分配不均，主要集中在5—9月，具有春秋季旱、夏季澇的特點[12]。研究區(qū)內(nèi)大部分區(qū)域為農(nóng)田，主要種植小麥、玉米等作物，核心建成區(qū)主要分布著蘆頭工業(yè)區(qū)、住宅區(qū)、商業(yè)區(qū)等。研究區(qū)內(nèi)煤炭、金、高嶺土、砂礦等礦產(chǎn)資源豐富，煤炭、采礦、鑄造和電鍍等工業(yè)發(fā)達，其中工業(yè)殘渣的堆放、廢液排放都會對周圍的環(huán)境產(chǎn)生一定的影響。由于當?shù)亟邓y以滿足農(nóng)業(yè)需求，經(jīng)污水處理廠處理后的生活廢水、工業(yè)污水經(jīng)由黃水河及永汶河后進入各支渠用于農(nóng)田灌溉，加劇農(nóng)田土壤中各類重金屬元素的富集，區(qū)域內(nèi)主要污染重金屬為Cr、Ni、Pb、Zn、Cu、Co、As、Cd。

1.2 樣品采集和理化分析

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

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

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

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

1.4 研究方法

1.4.1 富集系數(shù)

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

(1)

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

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

1.4.2 地理探測器模型

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

表1 交互探測表達式Table 1 Expression of interactive detection

2 結(jié)果與討論

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

研究區(qū)土壤樣品pH為7.00～7.85，平均值為7.43，屬于中性偏堿性土壤，土壤樣品中重金屬質(zhì)量濃度與富集系數(shù)統(tǒng)計見表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ū)土壤中重金屬受人為干擾較大，集聚趨勢比較明顯。本研究選取Mn作為參考元素，計算得到研究區(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，可以被認為是受到人為因素影響的主要重金屬元素，因此在研究區(qū)中富集現(xiàn)象明顯；Co、Cr、Ni富集系數(shù)低于1.5，且在研究區(qū)范圍內(nèi)標準偏差均小于等于0.2，說明這3種重金屬的富集主要受自然因素的影響較大，受到人為因素影響較小。

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

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

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

2.3 地理探測器模型分析

2.3.1 因子分析

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

圖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為例考察影響因子間對重金屬污染的交互作用。表4為影響因子間交互作用對Pb污染的解釋力，通過比較單一因子解釋力和交互作用解釋力的大小關(guān)系，判斷兩個影響因子的交互作用效果。經(jīng)比較，NDVI、SOM、DI、DR、DP、DEM、PMT間的交互作用均為非線性增強，也即影響因子間單因子解釋力合計小于影響因子交互作用的解釋力，Pb分布受到各影響因子間交互作用的影響。其中，DI和DP交互作用對Pb的解釋力達0.921，是影響Pb分布的主要因子，其他影響因子也會加強工業(yè)以及人口因素對Pb分布的影響。

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

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

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

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

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

圖3 地理探測器模型與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) 地理探測器因子探測結(jié)果表明，7種影響因子中DI對8種重金屬的平均解釋力最大,其次為SOM和PMT。不同重金屬的主要影響因子也存在不同，Cd、Pb主要歸因于DI，Zn、Cu主要歸因于SOM，As、Co、Ni主要歸因于PMT，Cr分布受到PMT、DP、SOM、DI的共同作用影響。

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

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