楊楠，連雅，，夏新，張明順，李宗超*

（1.中國環(huán)境監(jiān)測總站，北京 100012；2.北京建筑大學環(huán)境與能源工程學院，北京 100044）

土壤形成過程中受到母質(zhì)、地形、水文和農(nóng)業(yè)生產(chǎn)等自然因素和人類活動因素影響，在空間和時間尺度上均具有高度的異質(zhì)性[1]。在土壤污染調(diào)查中，為提高調(diào)查結(jié)果的精度，常加密布設點位[2]。但在后續(xù)跟蹤監(jiān)測中難以對全部點位開展長期監(jiān)測，因此需對地塊內(nèi)的調(diào)查點位進行篩選。WANG 等[3]、GAO 等[4]研究中和《地理信息 空間抽樣與統(tǒng)計推斷》（GB/Z 33451—2016）標準中提及可在均質(zhì)性區(qū)域內(nèi)篩選代表性點位。為獲得土壤重金屬均質(zhì)性區(qū)域，本文開展了相關(guān)分區(qū)方法研究，使地塊分區(qū)后子區(qū)域內(nèi)土壤重金屬含量趨于均質(zhì)化，子區(qū)域間趨于異質(zhì)化。

目前多數(shù)學者使用地統(tǒng)計學與GIS 結(jié)合的方法對土壤重金屬空間變異特征進行研究，利用半變異函數(shù)描述土壤重金屬空間變異結(jié)構(gòu)，使用基于半變異函數(shù)的克里金插值對污染空間進行預測，直觀地觀測重金屬的空間連續(xù)分布情況[5?9]。模糊聚類分區(qū)方法可同時考慮多個土壤屬性變量，進行土壤空間管理區(qū)域的劃分[10?13]，但該方法未考慮變量之間的相關(guān)性對聚類分析帶來的誤差影響[14]。主成分?模糊聚類分析方法，利用主成分分析提取原始數(shù)據(jù)的絕大部分信息并將多個具有相關(guān)性的指標轉(zhuǎn)化為少量互不相關(guān)的綜合性指標，再將主成分得分作為新指標進行聚類分析，極大提高了分區(qū)的效率和準確性，該方法目前被廣泛應用在基于有機質(zhì)等土壤養(yǎng)分指標的土壤空間均質(zhì)性分區(qū)研究中[15?17]，但將其應用于以重金屬含量為指標的土壤污染空間均質(zhì)性分區(qū)研究相對較少。

本研究以土壤重金屬均質(zhì)性分區(qū)為導向，利用地統(tǒng)計學與GIS 結(jié)合的方法對土壤重金屬元素的空間變異結(jié)構(gòu)特征進行描述，獲得重金屬空間插值結(jié)果；利用主成分?模糊聚類分析方法開展重金屬空間均質(zhì)性分區(qū)研究；通過方差分析驗證分區(qū)結(jié)果合理性，并通過與重金屬空間插值結(jié)果比較，討論分區(qū)結(jié)果可靠性。本研究旨在證明主成分?模糊聚類分析方法可用于重金屬均質(zhì)性分區(qū)，并為在均質(zhì)性區(qū)域內(nèi)篩選代表性點位奠定基礎(chǔ)。

1 材料與方法

1.1 研究區(qū)域概況

為充分研究主成分?模糊聚類分析分區(qū)方法的適用性，本文選擇了受自然和人為因素影響、土壤重金屬空間變異特征較復雜的某山區(qū)金屬礦下游的耕地地塊作為研究區(qū)域。該區(qū)域位于北回歸線以北，平均海拔為805 m，區(qū)域面積為5.54 km2，東西橫跨2 km，南北縱跨4 km。研究區(qū)域土壤由兩種成土母質(zhì)——泥盆系碎屑巖和泥盆系碳酸鹽巖發(fā)育而來，區(qū)域范圍內(nèi)土壤類型有紅壤和水稻土兩種，所占面積分別為3.62 km2和1.92 km2。

1.2 數(shù)據(jù)來源

將研究區(qū)域劃分為200 m×200 m 的網(wǎng)格，并在耕地內(nèi)布設樣點（112 個），如圖1 所示。通過GPS 開展野外定點采樣工作，采用單點采樣法，采樣深度為0～20 cm。采樣過程中避免與金屬工具接觸，樣品置于聚乙烯塑料袋中。土壤樣品經(jīng)風干、逐級研磨過100目篩后，開展8種重金屬檢測。采用鹽酸?硝酸?氫氟酸?高氯酸全消解的方法，Pb 和Cd 含量采用石墨爐原子吸收光譜法（GB 17141—1997）測定，Cr、Zn、Cu和Ni含量采用火焰原子吸收法（HJ 491—2009）測定，As 和Hg 含量采用原子熒光光譜法（GB 22105.1—2008和GB 22105.2—2008）測定。分析過程中均加入國家土壤標準物質(zhì)（GSS?6），重金屬回收率均在90%以上，所有樣品均設置3 個平行樣，平行樣的相對誤差小于10%，同時進行空白實驗，選10%的樣品做重復測試，相對誤差在±5%以內(nèi)。

1.3 研究方法

1.3.1 土壤重金屬含量地統(tǒng)計分析

基于地統(tǒng)計學對土壤重金屬含量的空間結(jié)構(gòu)特征進行探索性分析，包括空間自相關(guān)性分析和空間變異性分析。

半變異函數(shù)是地統(tǒng)計學函數(shù)，一般用來描述區(qū)域化變量的空間變異結(jié)構(gòu)[18]。數(shù)學計算公式為：

式中：N（h）為具有相同間隔距離h時的樣點對數(shù)量，Z（xi）和Z（xi+h）分別為位置在xi和xi+h上的變量實測值。半變異函數(shù)有4 個重要的參數(shù)：塊金、基臺值、塊金效應和變程。塊金值代表隨機變異?；_值代表由結(jié)構(gòu)變異和隨機變異組成的總體變異。塊金效應是塊金值與基臺值的比值，即隨機（人為）因素引起的空間變異占總體變異的比例，取值分別為：0%～25%、25%～50%、50%～75%、75%～100%時，代表的空間相關(guān)性程度分別為：強、明顯、中等、弱[19]。變程表示空間自相關(guān)的作用范圍，在變程范圍內(nèi)，變量具有空間自相關(guān)性，反之則不存在[20]。

克里金插值法是以變異函數(shù)理論和結(jié)構(gòu)分析為基礎(chǔ)，對區(qū)域化變量進行無偏最優(yōu)估計的插值方法。計算公式為：

式中：x0為待預測樣點，Z（xi）為樣點xi處重金屬的測試值，λi為樣點xi處測試值的權(quán)重。

1.3.2 主成分?模糊聚類分析

主成分分析可以通過正交變換將多個相關(guān)性變量轉(zhuǎn)化為幾個互不相關(guān)的綜合變量（主成分），一般提取特征值>1 的主成分進入后續(xù)數(shù)據(jù)分析與研究[21?22]。樣點主成分值（即主成分得分）的計算公式為：

式中：Fnj為第n個樣點第j項主成分的取值；m為常數(shù)，表示原始變量的個數(shù)；eji為第j項主成分第i項變量的載荷；fj為第j項主成分的特征值；Xni為第n個土壤樣點第i項變量的標準化值。

模糊聚類分析是根據(jù)樣本數(shù)據(jù)本身所具有的特征，采用相似性原則進行分類，最終使同類分區(qū)內(nèi)數(shù)據(jù)變異最小，而不同類分區(qū)間數(shù)據(jù)變異最大[23]。分類效果的優(yōu)劣主要由兩個指標[模糊效果指數(shù)（Fuzzy performancel index，F(xiàn)PI）和歸一化分類熵（Normalized classification entropy，NCE）]確定，計算公式分別為：

式中：c為分類數(shù)，n為樣本數(shù)，uik為第k個樣本屬于第i個聚類中心的模糊隸屬度。FPI和NCE取值均為0～1，F(xiàn)PI越接近0，聚類所得分區(qū)之間共用的數(shù)據(jù)越少，聚類效果就越好。NCE越接近0，分區(qū)內(nèi)的數(shù)據(jù)相似程度越高，聚類效果就越好。因此，F(xiàn)PI和NCE同時取最小值的聚類數(shù)為最佳分區(qū)數(shù)[24]。

1.3.3 數(shù)據(jù)分析和軟件制圖

使用SPSS 24.0 軟件進行描述性統(tǒng)計分析和主成分分析，使用GS+9.0 軟件進行半方差函數(shù)擬合，使用ArcGIS 10.2 軟件繪制重金屬含量空間分布預測圖和主成分得分空間分布預測圖，使用MZA 1.0.1 軟件進行模糊聚類分析，通過ArcGIS 10.2 軟件繪制分區(qū)結(jié)果圖。

2 結(jié)果與分析

2.1 土壤重金屬含量狀況

由柯爾莫戈洛夫?斯米諾夫（Kolmogorov?Smirnov）正態(tài)分布檢驗（P≤0.05）結(jié)果可知，Cr、Ni 含量符合正態(tài)分布，Cu對數(shù)變換后符合正態(tài)分布，Hg倒數(shù)變換后符合正態(tài)分布，As、Pb、Zn、Cd不符合正態(tài)分布。研究區(qū)域土壤重金屬含量統(tǒng)計情況見表1，由表1 可知，As、Cr、Pb、Zn、Cd、Cu、Ni、Hg 峰度均大于0，說明其數(shù)據(jù)分布曲線較為陡峭，為高狹峰。Cr 元素數(shù)據(jù)分布具有負偏性，說明Cr 元素含量中有偏低值；其余7 種元素數(shù)據(jù)分布具有正偏性，說明該7 種元素含量中均有偏高值，尤其除Ni 外，其他6 種元素偏度較高，說明其偏高的測試值較大。8 種重金屬變異程度由強到弱為：As>Pb>Zn>Cd>Hg>Cu>Ni>Cr，變異系數(shù)取值0%～10%為輕度變異，10%～100%為中度變異，大于100%為高度變異[25]，因此，Hg、Cu、Ni、Cr 為中度變異，As、Pb、Zn、Cd為高度變異。

表1 土壤重金屬含量統(tǒng)計特征描述Table 1 Description of statistical characteristics of heavy metal contents in soil

2.2 土壤重金屬含量地統(tǒng)計分析

2.2.1 土壤重金屬空間變異分析

地統(tǒng)計學一般要求所研究的區(qū)域化變量均服從正態(tài)或近似于正態(tài)分布[26]。本研究根據(jù)拉依達準則將樣本平均值加減3 倍標準差區(qū)間外的數(shù)據(jù)定為異常值，以平均值與3 倍標準差的和或差代替異常值[27]。對異常值處理后的各重金屬變量數(shù)據(jù)再次進行柯爾莫戈洛夫?斯米諾夫正態(tài)分布檢驗，Cr、Cu、Ni、Hg 4 個變量的數(shù)據(jù)符合正態(tài)分布，As、Pb、Zn、Cd 4 個變量的數(shù)據(jù)在進行對數(shù)變換后近似符合正態(tài)分布，均符合要求。

對研究區(qū)域As、Cr、Pb、Zn、Cd、Cu、Ni、Hg 8 種重金屬開展空間變異分析，對各重金屬數(shù)據(jù)構(gòu)建半變異函數(shù)模型，選擇決定系數(shù)較高且殘差最小的擬合模型作為最優(yōu)模型，模擬結(jié)果見表2。從擬合結(jié)果來看，8 種元素的決定系數(shù)均高于0.8，說明各變量的模型擬合精度較好，滿足后續(xù)分析要求。Zn、Cd元素的塊金效應低于25%，說明其具有強烈的空間自相關(guān)性；As、Cr、Pb、Cu、Ni、Hg 元素的塊金效應為25%～75%，說明其具有中等程度的空間自相關(guān)性。

表2 土壤重金屬變量半變異函數(shù)擬合及相關(guān)參數(shù)Table 2 Semivariogram fitting and related parameters of heavy metal variables in soil

8 種重金屬變程從大到小為：Cr>Pb>Hg>As>Cu>Zn>Ni>Cd，自相關(guān)范圍逐漸減小，除Cd 外（其有效變程為153 m），其余7種重金屬的變程均大于布設點位的影響范圍（200 m），說明土壤樣點基本在各金屬元素的空間自相關(guān)范圍內(nèi)，各元素的半變異函數(shù)模型能反應研究區(qū)域內(nèi)重金屬含量的空間變異情況。

2.2.2 土壤重金屬含量空間分布特征

根據(jù)各元素分布特征對數(shù)據(jù)進行相應的對數(shù)變換，根據(jù)趨勢分析結(jié)果剔除相應的全局趨勢，根據(jù)上述最佳半變異函數(shù)模型及其擬合參數(shù)，確定相應的插值參數(shù)，對8 種重金屬含量進行普通克里金插值，利用均方根誤差、標準均方根誤差和平均標準誤差作為指標，通過交叉驗證評估克里金插值模型的精度，確定最優(yōu)插值精度的插值模型，最終獲得各重金屬元素含量的空間分布預測圖，見圖2。

由圖2 可知，As、Pb、Zn、Cd、Cu、Hg 6 種重金屬含量的空間分布在整體上具有一定的相似性，大面積高值區(qū)在西部靠北區(qū)域，向東重金屬含量逐漸減??；Cr元素含量空間分布則與之相反，高值區(qū)在東部靠北區(qū)域，向西重金屬含量逐漸減小；而Ni 元素含量空間分布則無明顯的大范圍高值區(qū)域，只在西部靠北區(qū)域有小面積高值區(qū)。

2.3 土壤重金屬多變量均質(zhì)性分區(qū)研究

2.3.1 重金屬變量相關(guān)性分析

為使均質(zhì)性分區(qū)更符合研究區(qū)域?qū)嶋H重金屬含量分布情況，以重金屬含量原始測試值為數(shù)據(jù)源，進行主成分?模糊聚類分析和相關(guān)的分區(qū)結(jié)果討論。主成分分析前對指標間相關(guān)性進行分析，形成相關(guān)系數(shù)矩陣，見表3。8種重金屬中，除Ni與As、Cr、Pb之間相關(guān)系數(shù)未達顯著水平，其他指標間的相關(guān)性均達顯著水平（P<0.01），其中Cd 與Zn 的相關(guān)系數(shù)最大，達到0.989。Bartlett 球狀檢驗顯示Bartlett 值=1 319.933，P<0.05，指標間存在相關(guān)性，適合進行主成分分析；KMO檢驗顯示KMO值=0.784（>0.5），主成分分析的結(jié)果可接受。

表3 土壤重金屬變量相關(guān)系數(shù)矩陣Table 3 Correlation matrix of heavy metal variables in soil

2.3.2 主成分分析

對8種重金屬變量進行主成分分析，結(jié)果見表4。采用最大方差旋轉(zhuǎn)法可提取出3 個主成分——第一主成分（PC1）、第二主成分（PC2）和第三主成分（PC3），作為模糊聚類分析的3個指標，旋轉(zhuǎn)后其特征值分別為4.65（>1）、1.28（>1）和1.17（>1），方差貢獻率分別為58.13%、15.96%和14.62%，累計貢獻率為88.71%，即這3 個主成分涵蓋了原始數(shù)據(jù)信息總量的88.71%。根據(jù)公式（3）和表4中各成分載荷，計算112個樣點的第一、第二和第三主成分得分，在比較普通克里金法與反距離權(quán)重法交叉驗證精度后，選擇使用精度較高的普通克里金插值法對其進行插值，結(jié)果見圖3。

由表4 可知，第一主成分中，As、Pb、Zn、Cd、Cu、Hg 元素載荷較高，說明第一主成分主要反映了該區(qū)域As、Pb、Zn、Cd、Cu、Hg 的分布狀況，即第一主成分得分的空間分布與As、Pb、Zn、Cd、Cu、Hg 含量的空間分布呈現(xiàn)高度正相關(guān)。第二主成分中，Cr 元素載荷較高，說明第二主成分主要反映了Cr 的分布狀況，即第二主成分得分的空間分布與Cr 含量的空間分布呈現(xiàn)高度負相關(guān)。第三主成分中，Ni 元素載荷比較高，說明第三主成分主要反映了Ni 的分布狀況，即第三主成分得分的空間分布與Ni 含量的空間分布呈現(xiàn)高度正相關(guān)。

表4 土壤重金屬主成分分析結(jié)果Table 4 Principal component analysis results of soil heavy metals

2.3.3 聚類分析

以30 m×30 m 為柵格像元大小，從2.3.2 第一、第二和第三主成分得分的克里金插值圖中，將研究區(qū)域劃分為6 184 個柵格，分別提取每個柵格PC1、PC2 和PC3 得分3 個指標的數(shù)據(jù)，進行模糊聚類分析。軟件參數(shù)設置為[13，28]：最大迭代次數(shù)為300，收斂閾值為0.000 1，模糊指數(shù)為1.30，最大分區(qū)數(shù)為6，最小分區(qū)數(shù)為2，獲得NCE和FPI曲線如圖4 所示。由圖4 可知，將研究區(qū)域分為3 個分區(qū)是最佳的。分區(qū)結(jié)果如圖5所示。

2.3.4 分區(qū)結(jié)果合理性驗證

如圖5 所示，通過主成分?模糊聚類分析方法將研究區(qū)域分為3 個互相鑲嵌的子區(qū)域：分區(qū)1、分區(qū)2和分區(qū)3，土地面積占比分別為1%、30%和69%，所含樣點分別為2、33個和77個。為驗證分區(qū)結(jié)果的合理性，統(tǒng)計3 個分區(qū)各重金屬含量平均值并進行方差分析，結(jié)果見表5。方差分析結(jié)果顯示，8 個重金屬指標在3 個分區(qū)的均值均有顯著性差異（P<0.05）。3 個分區(qū)中，As、Pb、Zn、Cd、Cu、Ni、Hg 平均含量為分區(qū)1>分區(qū)2>分區(qū)3，Cr 平均含量為分區(qū)3>分區(qū)2>分區(qū)1，與土壤重金屬含量空間分布預測（圖2）所顯示的研究區(qū)域重金屬的空間分布基本一致。為進一步驗證分區(qū)結(jié)果，分別計算分區(qū)前后各重金屬的變異系數(shù)并進行比較（圖6），3 個分區(qū)內(nèi)8 個重金屬變量的變異系數(shù)，除Cr的分區(qū)1變異系數(shù)稍大于分區(qū)之前區(qū)域整體的變異系數(shù)，其他變量的分區(qū)變異系數(shù)均小于分區(qū)前整個研究區(qū)域。即分區(qū)后子區(qū)域內(nèi)的變異變小，而子區(qū)域之間的重金屬含量變異變大，分區(qū)合理。

表5 不同土壤重金屬分區(qū)平均含量統(tǒng)計及方差分析Table 5 Average content statistics and variance analysis of different soil heavy metal zones

3 討論

3 個分區(qū)互為鑲嵌的碎片化小區(qū)域在重金屬含量空間分布預測圖（圖2）中未得到全面展現(xiàn)，為對其進行具體討論，對碎片化小區(qū)域內(nèi)點位進行編號（圖5），分區(qū)1的碎片化區(qū)域內(nèi)點位為9號和12號，分區(qū)2的碎片化區(qū)域內(nèi)點位為1、2、3、4、5、6號和13號，分區(qū)3的碎片化區(qū)域內(nèi)點位為7、8、10號和11號，利用Anse?lin Local Moran′ I 指數(shù)對碎片化小區(qū)域的13 個點位的聚類及異常值情況進行分析。黑色9 號和12 號點位為As、Pb、Cd、Hg 含量影響的“高值聚類點位”，是研究區(qū)域重金屬含量最高值和次高值點位。紅色1、2、3、4、5、6 號和13 號點位為被As、Cd、Cu、Ni 含量影響的“被低值圍繞的高值異常點位”，所以其被分為分區(qū)2。而藍色7、8、10號和11號點位為As、Pb、Zn、Cu、Hg 含量影響的“高值聚類點位”，但是與其周圍的分區(qū)1 內(nèi)點位數(shù)據(jù)對比來說，屬于相對低值點位，所以被分為分區(qū)3。主成分?模糊聚類分析方法得到的分區(qū)與該研究區(qū)域?qū)嶋H土壤重金屬含量空間分布情況相符。而圖2 重金屬含量空間分布預測圖中未能具體展現(xiàn)這些點位的數(shù)據(jù)特征，可能是在地統(tǒng)計分析時由于異常值處理和（或）克里金插值的平滑效應被掩蓋了。碎片化小區(qū)域與其實際的大面積分區(qū)區(qū)域在空間上呈非連續(xù)分布，且通過上述討論，其內(nèi)的13 個點位重金屬含量呈現(xiàn)不同程度的“異?！保茄芯繀^(qū)域內(nèi)估計和預測重金屬含量的關(guān)鍵性點位，在進一步的點位優(yōu)化篩選研究前，可考慮優(yōu)先保留這13 個點位，同時將對應的碎片化區(qū)域設置為非優(yōu)化區(qū)域，除此之外的大面積分區(qū)區(qū)域作為最終分區(qū)進行研究。

通過主成分?模糊聚類分析方法得到的研究區(qū)域空間分區(qū)結(jié)果并不唯一，分區(qū)過程中存在不確定性因素。一是主成分的提取個數(shù)的不確定，主成分分析中主成分提取通?？紤]特征值（是否大于1）和累計貢獻率兩個因素，對于土壤重金屬污染的主成分分析，紀冬麗等[29]、劉慧琳等[30]采用最大方差旋轉(zhuǎn)法提取特征值大于1 的主成分進行分析，累計貢獻率均達85%以上，本研究中主成分提取方法和結(jié)果與之一致。二是主成分得分空間分布預測的不確定，對于土壤重金屬插值分析，肖艷桐等[31]指出應根據(jù)不同的數(shù)據(jù)分布和影響因素選擇適合的空間插值方法，本研究考慮到主成分得分呈對數(shù)正態(tài)分布，且空間自相關(guān)性強，同時為保證分區(qū)過程中的精度，通過比較普通克里金插值和反距離權(quán)重插值兩種方法的交叉驗證精度，最終選用精度相對較高的克里金插值對3 個主成分得分的空間分布進行預測。此外，李凱等[32]、謝云峰等[33]指出克里金插值易丟失局部細節(jié)，且克里金插值過程中，從數(shù)據(jù)檢驗、全局趨勢剔除到半變異函數(shù)擬合等，每一步的參數(shù)確定都會影響克里金插值結(jié)果，從而影響均質(zhì)性分區(qū)結(jié)果。因此分區(qū)過程中對數(shù)據(jù)進行準確分析并確定較為合適的參數(shù)對獲得可靠合理的分區(qū)結(jié)果來說至關(guān)重要。

4 結(jié)論

（1）Zn、Cd 含量具有強烈的空間自相關(guān)性，As、Cr、Pb、Cu、Ni、Hg 含量具有中等程度的空間自相關(guān)性，這為通過普通克里金插值方法獲得土壤重金屬空間分布預測圖提供了依據(jù)。

（2）利用主成分?模糊聚類分析方法將研究區(qū)域劃分為3 個均質(zhì)性分區(qū)，通過空間插值圖比較和方差分析對分區(qū)結(jié)果的合理性進行了驗證，除Cr 外，其他7 種元素均質(zhì)性變異系數(shù)均小于分區(qū)前整個研究區(qū)域，分區(qū)均合理。主成分?模糊聚類分析方法可對土壤重金屬空間變異特征較復雜的區(qū)域進行多變量均質(zhì)性分區(qū)，分區(qū)后不同分區(qū)之間重金屬含量差異顯著，分區(qū)內(nèi)各重金屬含量較為均勻，對區(qū)域土壤重金屬含量的準確估計與預測具有重要意義，同時為利用均質(zhì)性分區(qū)篩選代表性點位奠定基礎(chǔ)。