農(nóng)用地土壤污染調(diào)查分階段分區(qū)采樣布點*

許友靜， 張 清, 劉興旺, 李 峰, 葛 飛

(湘潭大學 環(huán)境與資源學院，湖南 湘潭 411105)

0 引言

區(qū)域土壤污染調(diào)查是風險評估、修復治理和科學利用的前提，而污染調(diào)查能否有效掌握污染物的空間分布信息，其采樣布點方法是關(guān)鍵[1-2].采樣往往會結(jié)合模型模擬預(yù)測與制圖來反映空間信息，因此采樣數(shù)量決定空間預(yù)測精度[3]，也會影響變量的模型預(yù)測結(jié)果[4]，采樣數(shù)量越多會產(chǎn)生越高的預(yù)測精度[5-6]，故采樣數(shù)量與預(yù)測精度的定量關(guān)系是優(yōu)化采樣策略的重要依據(jù)，對區(qū)域土壤污染調(diào)查具有重要的現(xiàn)實意義.

采樣布點研究是利用采樣點位結(jié)合變異函數(shù)和空間插值方法對未采樣區(qū)域進行預(yù)測，并根據(jù)預(yù)測誤差進行采樣數(shù)量的優(yōu)選[7].地統(tǒng)計條件模擬方法是基于概率分布的思路，同時獲得多個模擬結(jié)果，通過超過相關(guān)標準的概率閾值來確定污染區(qū)范圍[8].序貫高斯模擬(Sequential Gaussian Simulation, SGS)適用于性質(zhì)非均勻和不確定性的定量表征[9]，避免了克里格法存在的明顯平滑效應(yīng)[10-11].地統(tǒng)計學方法早已被廣泛應(yīng)用于土壤屬性的空間分布預(yù)測，用有限的樣本量對土壤屬性的空間分布進行數(shù)字制圖[12]，其中研究最多的是農(nóng)用地變量施肥，以土壤養(yǎng)分的空間變異為出發(fā)點和依據(jù)，根據(jù)有限的采樣點信息估算農(nóng)田土壤養(yǎng)分分布[13].相較于土壤養(yǎng)分，土壤污染受自然因素和人為因素的共同作用，具有更高的空間變異性[14].目前農(nóng)用地土壤污染調(diào)查的采樣策略通常將整個研究區(qū)作為一個整體，再結(jié)合地統(tǒng)計插值結(jié)果及預(yù)測精度優(yōu)化采樣數(shù)，對土壤污染物空間分布的高度差異性考慮不夠；在復雜地形區(qū)域僅以地形要素劃分采樣單元進行探討[13]；其他采樣單元的劃分多以經(jīng)驗判斷為主.

本研究基于地統(tǒng)計模擬和網(wǎng)格復合采樣探索一種農(nóng)用地土壤污染調(diào)查分階段分區(qū)采樣布點方法，第一階段利用初始網(wǎng)格和序貫高斯模擬(SGS)不確定性分析，按照污染物空間分布差異程度進行分區(qū)；第二階段分區(qū)布點，確定各區(qū)的合理采樣數(shù).利用收集的某鄉(xiāng)鎮(zhèn)農(nóng)用地土壤污染數(shù)據(jù)進行模擬及驗證，旨在提高農(nóng)用地土壤污染調(diào)查采樣布點效率，節(jié)省采樣成本.

1 材料與方法

1.1 研究區(qū)概況

研究區(qū)選取南方某鄉(xiāng)鎮(zhèn)，區(qū)域內(nèi)有一條主干河流呈西向東穿境，形成復雜的水系.地貌類型為丘陵，平均海拔70 m，亞熱帶季風氣候，年平均氣溫16～18 ℃.區(qū)域內(nèi)農(nóng)作物主要為水稻，農(nóng)用地土壤面積約159.94 km2，區(qū)域礦產(chǎn)資源豐富，由于長期開采導致農(nóng)用地土壤污染嚴重.本文收集研究區(qū)內(nèi)農(nóng)用地土壤長期跟蹤調(diào)查和采樣歷史數(shù)據(jù)，所有采樣點分布如圖1所示.

圖1 研究區(qū)采樣點空間分布Fig.1 Spatial distribution of sampling points in study area

1.2 序貫高斯模擬

序貫高斯模擬(SGS)是一種連續(xù)的變量序貫?zāi)M算法，可以提供多個可能的模擬實現(xiàn)，根據(jù)多個模擬實現(xiàn)，在模擬值接近真實值的前提下再現(xiàn)變量的空間變化特征，能根據(jù)多個模擬實現(xiàn)的平均值預(yù)測變量空間分布[15].SGS多次模擬實現(xiàn)之間的波動性為預(yù)測結(jié)果的不確定性提供了一個定量化的途徑[16]，通過局部不確定性的定量分析進行分區(qū)，空間某具體位置xm處污染物濃度的不確定性(即局部不確定性)可定義為在該位置1 000次模擬實現(xiàn)結(jié)果中污染物濃度低于給定濃度閾值c的概率：

(1)

式中：1 000為模擬次數(shù)；n(xm)為某具體空間位置xm處所產(chǎn)生的1 000個污染物濃度模擬實現(xiàn)結(jié)果中高于給定濃度閾值c的次數(shù).濃度閾值選取《土壤環(huán)境質(zhì)量 農(nóng)用地土壤污染風險管控標準(試行)》(GB15618—2018)中相應(yīng)污染物的篩選值或管控值.

1.3 分階段分區(qū)采樣策略

本研究并不進行實際采樣，而是根據(jù)已有采樣數(shù)據(jù)集，利用空間抽樣的思路模擬抽樣并進行驗證，提出的分階段分區(qū)采樣布點方法包括以下步驟：

(a) 初始網(wǎng)格間距的設(shè)定：本研究分階段分區(qū)采樣布點基于初始網(wǎng)格，由于地統(tǒng)計學采樣方法依賴于空間相關(guān)性[17]，空間相關(guān)性越強，分區(qū)效率越高[18]，故需首先確定最低有效的初始網(wǎng)格樣點進行分區(qū).設(shè)置1 km×1 km、2 km×2 km、3 km×3 km、4 km×4 km、5 km×5 km等5個間距的規(guī)則網(wǎng)格抽樣，抽取靠近網(wǎng)格中心位置的點位，采用莫蘭指數(shù)(Moran′s Ⅰ)評價得出合理的初始網(wǎng)格間距[19].

(b) 序貫高斯模擬不確定性分區(qū)及效果驗證：將初始網(wǎng)格點位通過SGS得到1 000次模擬結(jié)果，依次得到超過《土壤環(huán)境質(zhì)量 農(nóng)用地土壤污染風險管控標準(試行)》(GB15618—2018)中管控值和篩選值的概率分布，選取概率閾值0.8，將研究區(qū)劃分為污染物低值區(qū)(Ⅰ)、過渡區(qū)(Ⅱ)和高值區(qū)(Ⅲ)，SGS空間分布預(yù)測模擬過程通過GS+9.0軟件完成.分區(qū)后若各區(qū)污染物濃度平均值分化、變異系數(shù)變化明顯則說明分區(qū)結(jié)果比較合理[20]，若不發(fā)生顯著的變化則分區(qū)結(jié)果不理想.

(c) 分區(qū)確定合理采樣數(shù)及預(yù)測精度分析：各區(qū)設(shè)置0.4 km×0.4 km、0.5 km×0.5 km、0.6 km×0.6 km、0.7 km×0.7 km、0.8 km×0.8 km、0.9 km×0.9 km、1.0 km×1.0 km、1.5 km×1.5 km、2.0 km×2.0 km、3.0 km×3.0 km等10個規(guī)則網(wǎng)格繼續(xù)模擬抽樣，利用SGS平均實現(xiàn)預(yù)測污染物空間分布.利用10%點位作為獨立點驗證預(yù)測精度，得到采樣數(shù)與預(yù)測精度的變化關(guān)系，預(yù)測精度驗證過程及制圖在ArcGIS 10.5軟件完成.為了對比不分區(qū)的情況，將研究區(qū)作為一個區(qū)(О)進行合理采樣數(shù)及預(yù)測精度分析.

為了驗證預(yù)測精度，選用均方根誤差(Root Mean Square Error, RMSE)進行評價，RMSE值越低，預(yù)測精度越高[21-22]：

(2)

式中：n為驗證集樣點數(shù)；oi為驗證點位的觀測值；pi為驗證點位上的預(yù)測值.

考慮到各分區(qū)間污染物濃度差異較大，存在尺度/規(guī)模效應(yīng)，使用標準化均方根誤差(NRMSE)來比較各單元間的相對預(yù)測精度，較小的NRMSE值意味著更高的預(yù)測精度，NRMSE通過式(3)計算[13,23].

NRMSE=(RMSE/mean)×100

，

(3)

式中：mean為各樣點集的污染物平均值.

2 結(jié)果與討論

2.1 研究區(qū)污染物特征和初始網(wǎng)格確定

對研究區(qū)數(shù)據(jù)集pH和污染物濃度進行描述性統(tǒng)計分析，如表1所示.

表1 研究區(qū)數(shù)據(jù)集描述性統(tǒng)計

由表1可知，污染物變異系數(shù)(V)為79.31%，根據(jù)Wilding(1985)[19]提出的標準，CV≤15%屬于弱變異，15%35%為高度變異，研究區(qū)污染物為高度變異.峰度和偏度系數(shù)表明研究區(qū)污染物總體樣本數(shù)據(jù)不符合正態(tài)分布，經(jīng)對數(shù)變換后服從正態(tài)分布，符合地統(tǒng)計學分析的要求[24].圖2為5組初始網(wǎng)格點位的空間相關(guān)性定量分析結(jié)果.

由圖2可以看出，莫蘭指數(shù)值在0.5～0.88范圍內(nèi)，顯著性檢驗指標Z得分在7～25范圍內(nèi)，說明5組網(wǎng)格樣點的污染物含量都呈現(xiàn)出正的全局空間自相關(guān)性[25].但當網(wǎng)格間距增大時，莫蘭指數(shù)值和顯著性檢驗指標Z得分顯著降低，說明網(wǎng)格間距越大，點位空間自相關(guān)性越弱，同時顯著性程度也越低.網(wǎng)格間距為1 km×1 km和2 km×2 km時，空間自相關(guān)性較強且顯著性程度也較高，網(wǎng)格間距為3 km×3 km時空間自相關(guān)性較強，顯著性程度顯著降低，超過3 km×3 km后，空間自相關(guān)性和顯著性程度均顯著降低.綜合考慮相關(guān)性、顯著性和采樣成本，本研究選用2 km×2 km作為初始網(wǎng)格進行下一步的分區(qū)，2 km×2 km初始網(wǎng)格對應(yīng)119個初始樣點，符合相關(guān)研究報道[26]中準確、穩(wěn)定的變異函數(shù)需要觀測樣點數(shù)量達到50～150個的要求.

圖2 不同網(wǎng)格間距莫蘭指數(shù)值Fig.2 Moran’s index of different initial grid space

2.2 序貫高斯模擬不確定性分區(qū)

利用序貫高斯模擬不確定性進行分區(qū)，將研究區(qū)劃分為污染物低值區(qū)(Ⅰ)、過渡區(qū)(Ⅱ)和高值區(qū)(Ⅲ)，劃分結(jié)果如圖3所示.

圖3 研究區(qū)農(nóng)用地分區(qū)結(jié)果Fig.3 Agriculture sub-division result of the study area

樣點的代表性高低對精度影響較大，獲取代表性較高的采樣點能顯著提高精度，而將區(qū)域劃分為變量性質(zhì)較為均勻的分區(qū)后便于獲取代表性較高的采樣點[27].表2為各分區(qū)樣點污染物含量統(tǒng)計結(jié)果，低值區(qū)(Ⅰ)、過渡區(qū)(Ⅱ)和高值區(qū)(Ⅲ)分別包含樣本點位853個、331個、87個，平均值由分區(qū)前的0.58 mg/kg明顯分化為0.36 mg/kg、0.8 mg/kg、1.88 mg/kg，表明分區(qū)后各區(qū)平均值差異顯著，變異系數(shù)也由分區(qū)前的79.31%分別降低至36.11%、35.0%和26.6%，表明單元內(nèi)變異顯著降低，污染物濃度較為平均，分區(qū)取得較好的效果.

表2 各分區(qū)樣點污染物描述性統(tǒng)計

2.3 各區(qū)合理采樣數(shù)及預(yù)測精度分析

采樣點獲得的變量空間信息和代表性會對預(yù)測精度產(chǎn)生顯著影響，但存在飽和點，達到飽和點后再增加樣點并不能提高預(yù)測精度[28-29].這可能是因為采樣點在空間上聚集，使空間信息達到一定的飽和程度，額外增加的樣點代表性降低，即使采樣點數(shù)量增加，也無法提供進一步的空間信息[30].本研究將各區(qū)飽和樣點位數(shù)作為合理采樣數(shù)能獲得最佳效益，保證預(yù)測精度的同時避免冗余點位造成浪費.

圖4(a)、(b)、(c)、(d)為各分區(qū)10組網(wǎng)格采樣數(shù)與污染物空間分布預(yù)測精度的變化趨勢和曲線擬合.由圖4可知，隨著網(wǎng)格密度的下降，各分區(qū)RMSE值呈現(xiàn)先基本保持平穩(wěn)后升高的趨勢，表明預(yù)測精度隨著采樣密度下降而下降，并且存在拐點可以作為合理采樣數(shù)，這與Long J等[13]的研究結(jié)論一致.各區(qū)拐點有所不同，低值區(qū)(Ⅰ)、過渡區(qū)(Ⅱ)和高值區(qū)(Ⅲ)和研究區(qū)(О)拐點分別出現(xiàn)在2.0 km×2.0 km、0.9 km×0.9 km、0.6 km×0.6 km、0.9 km×0.9 km網(wǎng)格，合理采樣數(shù)分別為77個、118個、56個、452個，采樣密度分別為0.7個/km2、2.9個/km2、4.6個/km2、2.8個/km2.

圖4 各分區(qū)預(yù)測精度(RMSE)Fig.4 Prediction accuracy of each unit (RMSE)

圖5為各區(qū)采樣密度與相對預(yù)測精度關(guān)系，低值區(qū)(Ⅰ)、過渡區(qū)(Ⅱ)和高值區(qū)(Ⅲ)和研究區(qū)(О)NRMSE值分別在28.97%～35.43%、18.45%～25.41%、11.75%～20.33%、38.05%～44.69%之間，隨著采樣密度的降低，相對預(yù)測精度分別下降了6.46%、6.96%、8.58%、6.64%，高值區(qū)下降幅度最大，過渡區(qū)和研究區(qū)次之，低值區(qū)最低.預(yù)測精度變化幅度反映了預(yù)測精度對采樣密度變化的敏感度，變化幅度越大對采樣密度變化的敏感性越高[13]，高值區(qū)預(yù)測精度對采樣密度變化的敏感性最高，過渡區(qū)和研究區(qū)次之，低值區(qū)最低，故在污染物濃度越高的區(qū)域增加采樣點數(shù)更能使預(yù)測精度顯著提高.

圖5 各分區(qū)相對預(yù)測精度(NRMSE)Fig.5 Relative prediction accuracy of each unit (NRMSE)

進一步比較圖5發(fā)現(xiàn)，高值區(qū)標準化均方根誤差NRMSE最低，其次是過渡區(qū)和低值區(qū)，研究區(qū)最高，各區(qū)最低NRMSE分別約為12%、18%、28%和39%，表明高值區(qū)相對預(yù)測精度最好，過渡區(qū)和低值區(qū)次之，研究區(qū)最差.本研究中高值區(qū)、過渡區(qū)、低值區(qū)、研究區(qū)變異系數(shù)分別約為25%、35%、36%和79%，即空間變異最小的高值區(qū)取得了最好的相對預(yù)測精度，空間變異最大的整個研究區(qū)相對預(yù)測精度最差.多項研究表明空間變異對預(yù)測精度有負面影響[29-30]，在空間變異相對較小的區(qū)域會產(chǎn)生更好的預(yù)測精度，如Long J[13]通過對復雜區(qū)域不同地形下農(nóng)田土壤有機質(zhì)采樣密度與預(yù)測精度分析發(fā)現(xiàn)，空間變異最大的平原臺地NRMSE值最高，其次是丘陵地形，空間變異最低的河谷盆地NRMSE值最低，預(yù)測精度表現(xiàn)最好，與本研究結(jié)論一致.

3 結(jié)論

1) 研究區(qū)總體樣本數(shù)據(jù)污染物變異系數(shù)CV為79.31%，屬高度變異，經(jīng)對數(shù)變換后服從正態(tài)分布，符合地統(tǒng)計學分析的要求；研究區(qū)合適的初始網(wǎng)格為2 km×2 km，此時數(shù)據(jù)基本滿足空間自相關(guān)要求.

2) 序貫高斯模擬不確定性分區(qū)后平均值由分區(qū)前的0.58 mg/kg明顯分化為低值區(qū)(Ⅰ) 0.36 mg/kg、過渡區(qū)(Ⅱ) 0.8 mg/kg、高值區(qū)(Ⅲ) 1.88 mg/kg，變異系數(shù)也由分區(qū)前的79.31%分別降低至低值區(qū)(Ⅰ)36.11%、過渡區(qū)(Ⅱ)35.0%和高值區(qū)(Ⅲ) 26.6%，分區(qū)方法取得較好的效果.

3) 分區(qū)確定合理采樣數(shù)為低值區(qū)(Ⅰ) 77個、過渡區(qū)(Ⅱ) 188個和高值區(qū)(Ⅲ) 56個，對應(yīng)NRMSE分別為28%、18%、12%，相比不進行分區(qū)(О)時合理采樣數(shù)452個和NRMSE 39%，以較少的采樣點取得了較好的預(yù)測精度.

4) 對比各區(qū)采樣密度和預(yù)測精度發(fā)現(xiàn)，污染物濃度越高的分區(qū)預(yù)測精度對采樣密度越敏感，污染物濃度越高的區(qū)域增加采樣點數(shù)更能使預(yù)測精度顯著提高；污染物空間變異越小的分區(qū)呈現(xiàn)越高的相對預(yù)測精度.本文提出的分階段分區(qū)采樣方法能顯著減少點位，節(jié)省采樣成本，提高采樣效率，為管理部門農(nóng)用地土壤污染調(diào)查采樣布點設(shè)計提供參考.