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

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

0 引言

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

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

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

1 材料與方法

1.1 研究區(qū)概況

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

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

1.2 序貫高斯模擬

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

(1)

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

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

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

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

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

(c) 分區(qū)確定合理采樣數(shù)及預(yù)測(cè)精度分析：各區(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個(gè)規(guī)則網(wǎng)格繼續(xù)模擬抽樣，利用SGS平均實(shí)現(xiàn)預(yù)測(cè)污染物空間分布.利用10%點(diǎn)位作為獨(dú)立點(diǎn)驗(yàn)證預(yù)測(cè)精度，得到采樣數(shù)與預(yù)測(cè)精度的變化關(guān)系，預(yù)測(cè)精度驗(yàn)證過程及制圖在ArcGIS 10.5軟件完成.為了對(duì)比不分區(qū)的情況，將研究區(qū)作為一個(gè)區(qū)(О)進(jìn)行合理采樣數(shù)及預(yù)測(cè)精度分析.

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

(2)

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

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

NRMSE=(RMSE/mean)×100

，

(3)

式中：mean為各樣點(diǎn)集的污染物平均值.

2 結(jié)果與討論

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

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

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

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

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

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

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

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

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

樣點(diǎn)的代表性高低對(duì)精度影響較大，獲取代表性較高的采樣點(diǎn)能顯著提高精度，而將區(qū)域劃分為變量性質(zhì)較為均勻的分區(qū)后便于獲取代表性較高的采樣點(diǎn)[27].表2為各分區(qū)樣點(diǎn)污染物含量統(tǒng)計(jì)結(jié)果，低值區(qū)(Ⅰ)、過渡區(qū)(Ⅱ)和高值區(qū)(Ⅲ)分別包含樣本點(diǎn)位853個(gè)、331個(gè)、87個(gè)，平均值由分區(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ū)樣點(diǎn)污染物描述性統(tǒng)計(jì)

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

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

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

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

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

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

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

3 結(jié)論

1) 研究區(qū)總體樣本數(shù)據(jù)污染物變異系數(shù)CV為79.31%，屬高度變異，經(jīng)對(duì)數(shù)變換后服從正態(tài)分布，符合地統(tǒng)計(jì)學(xué)分析的要求；研究區(qū)合適的初始網(wǎng)格為2 km×2 km，此時(shí)數(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個(gè)、過渡區(qū)(Ⅱ) 188個(gè)和高值區(qū)(Ⅲ) 56個(gè)，對(duì)應(yīng)NRMSE分別為28%、18%、12%，相比不進(jìn)行分區(qū)(О)時(shí)合理采樣數(shù)452個(gè)和NRMSE 39%，以較少的采樣點(diǎn)取得了較好的預(yù)測(cè)精度.

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