農(nóng)用地土壤重金屬污染調(diào)查最優(yōu)網(wǎng)格尺度及布點優(yōu)化方法

秦旭芝，鄭涵文，何瑞成，丁文娟，張軍,崔冠男，史沛麗，謝云峰*

1.廣西壯族自治區(qū)生態(tài)環(huán)境監(jiān)測中心 2.桂林理工大學(xué)環(huán)境科學(xué)與工程學(xué)院 3.生態(tài)環(huán)境部土壤與農(nóng)業(yè)農(nóng)村生態(tài)環(huán)境監(jiān)管技術(shù)中心

土壤重金屬污染來源多樣，由于重金屬具有難降解、難遷移等特點，一旦形成污染會在土壤中長期存在，并且在空間分布上呈現(xiàn)出明顯的異質(zhì)性和不連續(xù)性，加大準確判斷土壤中污染物空間分布的難度[1]?，F(xiàn)有技術(shù)條件下缺乏直接確定土壤污染物空間分布的方法，只能基于采樣調(diào)查數(shù)據(jù)通過數(shù)學(xué)模型進行空間插值[2-4]，而這種研究方法的準確性受前期采樣調(diào)查結(jié)果的影響很大。大量研究表明，對于空間分布變異性很大的污染物，布點方案是影響污染調(diào)查結(jié)果準確性的重要因素，如Jenkins等[5]對土壤中三硝基甲苯空間分布變異性的研究發(fā)現(xiàn)，95%的變異性來源于采樣位置；Argyraki等[6]發(fā)現(xiàn)由采樣位置造成的不確定性占整個污染物含量估計的不確定性的50%。結(jié)合調(diào)查區(qū)域的歷史背景、地理環(huán)境，科學(xué)有效地制定布點方案能大幅提高土壤污染調(diào)查的準確性[7-9]。目前常用的布點方法主要包括兩大類[10]：1)基于模型的布點方案，如網(wǎng)格布點法、隨機布點法等[11]；2)基于設(shè)計的布點方案，研究者基于調(diào)查區(qū)域的環(huán)境背景及土地利用歷史等對可能存在污染的區(qū)域進行加密布點[12]。Brus等[13-15]對2種布點方案進行了詳細論述和比較。

實際案例中多種布點方案常協(xié)同使用[16-20]，但如何確定最優(yōu)采樣密度是土壤污染調(diào)查的重點和難點問題[21-24]。謝志宜等[25]在開展珠三角區(qū)域耕地土壤質(zhì)量研究時發(fā)現(xiàn)，布點密度越大，調(diào)查區(qū)域土壤質(zhì)量刻畫越細致，但調(diào)查成本也越高，8 km×8 km的網(wǎng)格，能較好地保證調(diào)查結(jié)果的準確性，同時也能夠有效降低調(diào)查成本。工業(yè)場地土壤通常具有重金屬污染嚴重、污染重金屬單一且分布較為均勻等特點[26-29]，重金屬污染調(diào)查尺度較小，調(diào)查時通常先采用系統(tǒng)布點法，排查整個調(diào)查區(qū)域內(nèi)可能存在的土壤環(huán)境風(fēng)險，然后結(jié)合分區(qū)布點法和專業(yè)判斷法對可能存在環(huán)境風(fēng)險的區(qū)域進行重點調(diào)查。蘭鵬鵬等[30]以40 m×40 m網(wǎng)格對湖南某冶煉場地進行分區(qū)布點，結(jié)果發(fā)現(xiàn)場地污染主要在1.0～2.5 m深層空間內(nèi)。

為確保調(diào)查區(qū)域內(nèi)污染物空間分布的準確性，常需采用分階段的方式對研究區(qū)域進行調(diào)查，主要包括基礎(chǔ)信息調(diào)查、初步調(diào)查和加密調(diào)查等環(huán)節(jié)。場地污染物空間分布調(diào)查的誤差主要來源于污染物空間分布插值結(jié)果不確定性較大的區(qū)域[31-34]，污染程度高估會導(dǎo)致調(diào)查費用和修復(fù)費用的增加，污染程度低估會導(dǎo)致環(huán)境風(fēng)險無法得到有效控制。點位優(yōu)化的目的在于在尋求降低調(diào)查結(jié)果不確定性的最佳樣本量，同時有效降低調(diào)查成本。本研究旨在探索農(nóng)用地土壤重金屬污染調(diào)查最優(yōu)布點密度及布點優(yōu)化，以提高土壤污染調(diào)查方法對污染區(qū)范圍和污染程度的估計精度，并為土壤污染調(diào)查提供方法學(xué)支持。

1 材料與方法

1.1 研究區(qū)域

研究區(qū)域位于中南某老工業(yè)園區(qū)周邊的農(nóng)用地，地處111°58′E～113°05′E，27°20′55″N～28°05′40″N，調(diào)查區(qū)域面積約75萬m2，所在區(qū)域總體地貌輪廓北、西、南地勢高，中部、東部地勢低平，但地勢起伏較為和緩，高程為52.17～79.69 m，高差不大，地貌類型多樣，山地、丘陵、崗地、平原、水面俱備(圖1)。區(qū)域四季分明，降水充沛，盛夏高溫，冬季寒冷。年均降水量為 1 425 mm，4—7月降水較集中，期間多有洪水發(fā)生，年均氣溫為17.5 ℃。冬季多西北風(fēng)，夏季多東南風(fēng)。研究區(qū)域周邊的工業(yè)園區(qū)具有50多年的化工生產(chǎn)歷史，工業(yè)企業(yè)數(shù)量多，“三廢”排放量大，污染物種類多，對周邊農(nóng)用地土壤環(huán)境影響較大。

注： 紅色表示高程較高；黃色表示高程中等；綠色表示高程較低。圖1 研究區(qū)域內(nèi)地形地勢示意Fig.1 Map of terrain within the scope of study

1.2 土壤采集與分析

按照70 m×70 m網(wǎng)格設(shè)計要求，共布設(shè)168個表層土壤(0～20 cm)采樣點，現(xiàn)場利用GPS獲取采樣點經(jīng)緯度，使用木鏟進行表層土壤采樣，低溫保存于密封袋中帶回實驗室。土壤樣品經(jīng)風(fēng)干去除雜物后，研磨過100目尼龍篩。采用HNO3-H2O2消解法進行樣品預(yù)處理，采用石墨爐原子吸收光譜法測定土壤Cd濃度。每個樣品設(shè)3組平行樣，選10%的樣品做重復(fù)測試，并進行空白試驗，測定過程中添加國家標(biāo)準物質(zhì)(GSS-6)進行質(zhì)量控制。同時測定土壤理化性質(zhì)(表1)。另外還檢測了其他5種重金屬(Cu、Ni、Pb、Hg、Zn)濃度。檢測結(jié)果顯示，除Ni外其他重金屬均存在超標(biāo)情況，其中Cd超標(biāo)點位最多，其他重金屬僅個別點位超標(biāo)。

表1 土壤理化性質(zhì)參數(shù)Table 1 Physical and chemical properties of soil

1.3 子樣本的抽取

將調(diào)查區(qū)域分為70 m×70 m、100 m×100 m、160 m×160 m和200 m×200 m 4種網(wǎng)格尺度。先隨機抽取18個點位作為獨立驗證集〔圖2(a)〕，構(gòu)建除70 m×70 m尺度外的其他3種尺度網(wǎng)格，將距離網(wǎng)格中心點最近的點位作為該尺度下的采樣點，最終得到100 m×100 m尺度網(wǎng)格點位數(shù)87個，160 m×160 m 尺度網(wǎng)格點位數(shù)37個，200 m×200 m尺度網(wǎng)格點位數(shù)24個(圖2)。

圖2 抽取樣本子集及其驗證點位分布Fig.2 Sampling subset and distribution of verification points

1.4 土壤Cd濃度空間插值及其精度評價

利用ArcGIS 10.5軟件對不同調(diào)查尺度下土壤Cd濃度進行空間分析，利用Matheron矩估計法，以權(quán)重最小二次方擬合理論半方差函數(shù)模型。采用克里格(Kriging)法對4種采樣尺度測定結(jié)果進行插值，分析不同采樣尺度檢測結(jié)果的精確度。由于4種采樣尺度土壤中Cd濃度均不符合正態(tài)分布，使用BOX-COX變換進行數(shù)據(jù)轉(zhuǎn)換，經(jīng)K-S 檢驗后符合正態(tài)分布。評價結(jié)果采用半方差函數(shù)及交叉驗證確定最佳精確度。結(jié)合調(diào)查區(qū)域土壤pH，根據(jù)GB 15618—2018《土壤環(huán)境質(zhì)量 農(nóng)用地土壤污染風(fēng)險管控標(biāo)準(試行)》，采用0.6 mg/kg作為區(qū)域內(nèi)Cd污染評價標(biāo)準。

2 結(jié)果與評價

2.1 不同采樣尺度下土壤Cd濃度變異性

對4種采樣尺度下土壤Cd進行半變異分析,結(jié)果如圖3和表2所示。從圖3看出，4種采樣尺度的采樣點間距離均小于空間相關(guān)距離，說明4種采樣尺度均適用于空間變異分析，但不同尺度下的空間相關(guān)距離存在差異。從圖3還可看出，土壤Cd濃度均表現(xiàn)出了一定的線性趨勢，說明存在較好的半方差結(jié)構(gòu)，且采樣尺度越小線性趨勢越明顯。塊金值通常被認為是代表隨機變異的量，而基臺值代表變量空間變異的結(jié)構(gòu)性方差，二者比值為塊金系數(shù)。當(dāng)塊金系數(shù)小于25%時，表示變量具有強烈的空間相關(guān)性；塊金系數(shù)為25%～50%，表示變量具有明顯的空間自相關(guān)；塊金系數(shù)為50%～75%，表示變量具有中等的空間自相關(guān)；塊金系數(shù)大于75%時，表示變量空間自相關(guān)性較小。從表2可知，70 m×70 m、100 m×100 m、160 m×160 m、200 m×200 m 4種采樣尺度下土壤Cd濃度的塊金系數(shù)分別為65.82%、46.40%、58.65%和45.37%，均在25%～75%之間，表明研究區(qū)土壤Cd濃度具有中等空間相關(guān)性，受結(jié)構(gòu)性因素和隨機因素影響。

圖3 不同采樣尺度下土壤Cd濃度半方差分布Fig.3 Semi-variance distribution of Cd content in soils at different sampling scales

表2 不同采樣尺度下土壤Cd半變異函數(shù)參數(shù)Table 2 Semi-variant function parameters of Cd in soils at different sampling scales

2.2 土壤Cd濃度空間預(yù)測精度比較

將4種網(wǎng)格尺度下土壤Cd濃度克里格插值進度交叉驗證(Cross Validation)，結(jié)果見表3。由表3可見，總體而言4種采樣尺度下采樣網(wǎng)格越小，平均誤差(ME)絕對值越小。70 m×70 m和100 m×100 m 采樣尺度的均方根誤差(RMSE)和平均標(biāo)準誤差(ASE)較小且接近，表明70 m×70 m和100 m×100 m網(wǎng)格采樣精度接近，且能較好地滿足預(yù)測要求，預(yù)測值與真實值相近。

表3 不同采樣尺度下克里格法插值精度估計Table 3 Estimation of Kriging interpolation accuracy at different sampling scales

圖4為不同采樣網(wǎng)格尺度下克里格插值方法土壤Cd濃度的獨立驗證結(jié)果。由圖4可知，網(wǎng)格密度越小，克里格插值預(yù)測精度越高，且4種采樣密度下|RMSE-ASE|隨采樣網(wǎng)格尺度的增大而增大，說明采樣網(wǎng)格越小，插值預(yù)測的結(jié)果越接近真實值。200 m×200 m網(wǎng)格尺度下精度遠小于其他較小的網(wǎng)格尺度?？傮w而言，采樣網(wǎng)格越大，土壤Cd的預(yù)測越平滑，對細節(jié)反映能力越弱，結(jié)合半方差函數(shù)分析及ME、RMSE、ASE 3項指標(biāo)，100 m×100 m網(wǎng)格能較為準確地反映場地Cd濃度分布狀況。

圖4 不同采樣網(wǎng)格尺度下土壤Cd濃度克里格 插值的獨立驗證結(jié)果Fig.4 Independent validation results of Kriging interpolation of Cd content in soils at different sampling grid scales

2.3 污染范圍預(yù)測不確定性分析

污染范圍界定得準確與否以及對預(yù)測不確定性區(qū)域的量化對場地的修復(fù)治理邊界以及相關(guān)決策的制定具有重要影響，污染范圍界定過大，將實際未污染的區(qū)域插值預(yù)測為污染區(qū)域，會過多地耗費人力、物力，加大修復(fù)成本，造成一定的經(jīng)濟負擔(dān)；而污染范圍界定過小，將實際污染的區(qū)域插值預(yù)測為非污染區(qū)域，導(dǎo)致場地污染風(fēng)險不能有效降低或消除。

利用污染預(yù)測標(biāo)準誤差圖，量化不同空間位置預(yù)測結(jié)果的不確定性。結(jié)果顯示，采樣點周邊位置的誤差通常很小，預(yù)測標(biāo)準誤差大的區(qū)域主要集中在樣點較為稀疏且有高值點的區(qū)域。為進一步明確預(yù)測結(jié)果的不確定性區(qū)域，將污染范圍預(yù)測圖和預(yù)測標(biāo)準誤差圖進行空間疊加計算，假定將預(yù)測值與誤差值之差大于0.6 mg/kg的范圍劃分為污染區(qū)域，二者之和小于0.6 mg/kg的范圍劃為未污染區(qū)域，其他為不確定性區(qū)域，劃定結(jié)果見圖5。由圖5可見，不確定區(qū)域主要分布在西南側(cè)和北側(cè)區(qū)域，這2個區(qū)域從20世紀60年代至21世紀初先后存在多家從事化工生產(chǎn)活動的企業(yè)，因生產(chǎn)條件有限、排污過程的跑冒滴漏會導(dǎo)致地塊土壤受到污染，污染物在地下水水力作用的影響下遷移擴散，使得地塊周邊農(nóng)田土壤成片受到污染，且西南側(cè)和北側(cè)區(qū)域內(nèi)還存在池塘等使前期該區(qū)域采樣點較少，即不確定區(qū)域主要存在樣點稀疏且高值區(qū)域。

圖5 土壤Cd污染調(diào)查點位分布Fig.5 Distribution of samples for investigation of Cd pollution in soils

2.4 不確定性區(qū)域加密驗證

通過預(yù)測標(biāo)準誤差圖劃定不確定性區(qū)域，按初步調(diào)查超標(biāo)點位鄰近原則，在不確定區(qū)域內(nèi)的初步調(diào)查超標(biāo)點位周邊選取最近點位作為詳細調(diào)查點位(圖5)，對初步調(diào)查和詳細調(diào)查Cd濃度結(jié)果進行統(tǒng)計分析，結(jié)果見表4。從表4可知，100 m×100 m網(wǎng)格初步調(diào)查結(jié)果空間位置精確度為78.89%，對場地污染范圍預(yù)測不確定性影響較大的主要是高估區(qū)域，占比為14.23%，在不確定性區(qū)域加密布點以后，空間位置精確度提高到86.25%，提高了7.36個百分點，且污染高估區(qū)域占比明顯下降。對污染物濃度統(tǒng)計分析發(fā)現(xiàn)，100 m×100 m網(wǎng)格初步調(diào)查與詳細調(diào)查結(jié)果差距不大，與所有樣點進行對比發(fā)現(xiàn)，100 m×100 m網(wǎng)格Cd平均濃度略大于70 m×70 m網(wǎng)格，但二者誤差滿足HJ/T 166—2004《土壤環(huán)境監(jiān)測技術(shù)規(guī)范》所規(guī)定的25%范圍。圖6為土壤Cd污染范圍預(yù)測結(jié)果。由圖6可見，Cd污染主要集中在場地西側(cè)，100 m×100 m網(wǎng)格加密前和加密后的差異主要集中在不確定性區(qū)域范圍內(nèi)。100 m×100 m 網(wǎng)格精確度有明顯提升，同時調(diào)查點位數(shù)量下降了35%。

表4 土壤Cd統(tǒng)計特征Table 4 Statistical characteristics of Cd in soil heavy metals

圖6 土壤Cd污染范圍預(yù)測結(jié)果Fig.6 Prediction results of Cd pollution range in soils

本研究旨在探索科學(xué)合理的農(nóng)用地土壤污染調(diào)查最佳布點密度及點位優(yōu)化方法，在分析土壤污染物空間自相關(guān)性的基礎(chǔ)上，利用地統(tǒng)計學(xué)方法進行空間分布插值，基于預(yù)測誤差和污染評價標(biāo)準，劃定調(diào)查結(jié)果不確定性區(qū)域，結(jié)合初步調(diào)查超標(biāo)點位和地形、地貌等條件，在不確定性區(qū)域中加密布點，從而提高調(diào)查精度。該方法依據(jù)地統(tǒng)計學(xué)基本假設(shè)，適用于空間分布具有中等相關(guān)性及以上的污染物，胡克林[35-37]等研究表明，表現(xiàn)出較明顯空間相關(guān)性，但對于空間自相關(guān)性較弱，僅在局部區(qū)域個別點位超標(biāo)的重金屬、多環(huán)芳烴等污染物，該方法的適用性需進一步驗證，這與謝云峰等[1]的研究結(jié)果類似。本研究采樣布點方法與目前概率加密布點法和模擬退火算法的各類加密方法相似[38-39]，基于初步調(diào)查點位數(shù)據(jù)進行地統(tǒng)計學(xué)分析，判斷污染物是否具有空間自相關(guān)性，對于具有空間自相關(guān)性的污染地塊再進行參數(shù)化計算，從而優(yōu)化加密布點區(qū)域，做到既有效降低調(diào)查成本，又能提高調(diào)查精度。

3 結(jié)論

(1)對于采用基于模型化采樣布點的污染場地調(diào)查而言，通過獨立驗證和交叉驗證模型方法，能夠?qū)Σ蓸泳W(wǎng)格尺度進行有效優(yōu)化，本研究中100 m×100 m網(wǎng)格能夠在準確反映場地污染程度及范圍的前提下，降低采樣及分析成本。

(2)基于空間插值預(yù)測標(biāo)準誤差與預(yù)測值進行運算，能合理確定場地污染的不確定性范圍，通過對不確定性區(qū)域進行針對性加密布點，能顯著提高預(yù)測精確度。