安永龍，萬利勤，李 霞，殷志強(qiáng)，衛(wèi)曉峰，何澤新，賈鳳超

(1.中國(guó)地質(zhì)環(huán)境監(jiān)測(cè)院，北京 100081；2.北京礦產(chǎn)地質(zhì)研究院，北京 100012)

土壤是成土母巖經(jīng)過風(fēng)化、剝蝕、搬運(yùn)、沉積等地質(zhì)作用后疊加成土作用形成的自然物質(zhì)體。在成土的每個(gè)階段，都可能會(huì)受到外界環(huán)境(氣候、海拔、地貌、人類活動(dòng)等[1～2])的作用，因此在區(qū)域化土壤的物理化學(xué)性質(zhì)和空間序列中形成了特有的空間變異性。傳統(tǒng)的統(tǒng)計(jì)學(xué)方法往往忽略了土壤樣本的固有屬性數(shù)據(jù)和空間屬性數(shù)據(jù)之間的交互式影響[3]，導(dǎo)致元素的異質(zhì)性和空間依賴性無法很好地表達(dá)[4]。目前已有不少國(guó)內(nèi)外學(xué)者針對(duì)不同區(qū)域時(shí)空尺度，利用地統(tǒng)計(jì)學(xué)結(jié)合GIS的方法解決土壤元素空間變異性問題。如崔萌等[5]運(yùn)用地統(tǒng)計(jì)學(xué)和主成分分析法對(duì)北京市大桃主產(chǎn)區(qū)土壤5種重金屬含量和空間結(jié)構(gòu)特征進(jìn)行了研究，結(jié)果顯示Pb、As、Cr、Hg含量受結(jié)構(gòu)性因素影響，而Cd含量受結(jié)構(gòu)性因素和隨機(jī)性因素共同影響；汪璇等[6]運(yùn)用地統(tǒng)計(jì)學(xué)方法對(duì)三峽庫(kù)區(qū)表層土壤部分微量元素進(jìn)行了空間變異特征分析，結(jié)果表明7種微量元素為中等自相關(guān)性，2種為弱自相關(guān)性，且4個(gè)方向上變異程度不明顯；Roger等[7]通過環(huán)境因子作為數(shù)據(jù)對(duì)土壤養(yǎng)分指數(shù)完成克里格插值，表明克里格插值與線性回歸共同分析可以一定程度提高土壤養(yǎng)分的預(yù)測(cè)精度。

承德作為北京的“菜籃子”基地，土壤安全直接關(guān)系著首都人民的健康，在這樣嚴(yán)峻的背景下摸清該地區(qū)土壤重金屬本底及空間分布狀況具有重要意義[8-9]。已有部分學(xué)者在區(qū)內(nèi)按不同流域開展了土壤重金屬調(diào)查研究工作，如孫厚云等[10]對(duì)承德市灤河流域12種重金屬地球化學(xué)基線值進(jìn)行厘定，發(fā)現(xiàn)表層土壤Cd、V、Ti等基線值高于河北省背景值，且Cd、Pb、Cu的累積程度較大；王夢(mèng)雨等[11]對(duì)承德市柳河流域農(nóng)田土壤重金屬進(jìn)行調(diào)查，發(fā)現(xiàn)Cd、Zn、Ni存在不同程度積累現(xiàn)象。本次研究基于地統(tǒng)計(jì)學(xué)和GIS技術(shù)，對(duì)承德市全域范圍內(nèi)表層土壤重金屬As、Cd、Cr、Cu、Hg、Ni、Pb、Zn等的空間結(jié)構(gòu)及分布特征進(jìn)行解析，并對(duì)土壤重金屬與Se之間的相關(guān)性進(jìn)行了探索研究，旨在掌握承德地區(qū)表層土壤重金屬空間分布狀況及其特征，并進(jìn)一步判斷重金屬與Se空間上是否存在一定關(guān)聯(lián)，為承德地區(qū)著力打造健康安全的生態(tài)特色農(nóng)果產(chǎn)業(yè)鏈，支撐“承德樣板”[12]落地，開展國(guó)土空間規(guī)劃和用途管控工作提供基礎(chǔ)支撐。

1 研究區(qū)概況

承德市處于華北和東北兩個(gè)地區(qū)的過渡地帶，南鄰北京，東靠遼寧，北接內(nèi)蒙古高原，地理位置十分重要。地勢(shì)北高南低，北部為壩上高原地貌，中南部為山地、丘陵地貌，地貌展布具有獨(dú)特的地理優(yōu)勢(shì)，為特色農(nóng)果產(chǎn)品和道地藥材種植提供了有利的地域條件，杏仁、山楂、板栗、金蓮花、黃芩等都是承德的名優(yōu)特產(chǎn)。承德市河流分屬三個(gè)流域四大水系，即灤河流域、三河流域、遼河流域，年產(chǎn)水量37.6×108m3，是京津唐的重要供水源地。區(qū)內(nèi)屬溫帶大陸性季風(fēng)型半濕潤(rùn)山地氣候，四季分明，雨熱同期。年平均溫度5.6 ℃，無霜期60～180 d，年降水量451～850 mm。土壤類型以褐土、棕壤、潮土、灰色森林土、粗骨土、栗鈣土等為主。

承德礦產(chǎn)資源分布較多，以鐵礦資源最為豐富，是我國(guó)僅次于攀枝花的第二大釩鈦磁鐵礦資源基地[13]。鐵礦和金礦集中分布于雙灤區(qū)和灤平縣北部以及寬城滿族自治縣周邊，銀礦、鉬礦、銅礦等分布較為分散。

圖1 研究區(qū)土壤樣點(diǎn)分布Fig.1 Distribution of soil sampling sites in the study area

2 材料與方法

2.1 樣品采集與分析方法

本次研究利用了2017年承德地區(qū)土地質(zhì)量地球化學(xué)調(diào)查數(shù)據(jù)，采樣點(diǎn)分布情況見圖1。土壤樣品采集過程中采用GPS定位結(jié)合地形圖定點(diǎn)方式，充分考慮土壤類型、土地利用方式。以土地利用類型為采樣的基本單元，土壤采樣深度為0～20 cm，并以“S”或“N”形采樣，去除雜草、礫石、蟲殼、糞便等雜物，一個(gè)組合樣由3～5個(gè)子樣點(diǎn)等組分混合而成，子樣點(diǎn)與采樣單元內(nèi)采集的土壤類型一致。采集樣品的原始重量不低于1 kg，過20目篩后樣品重量不少于500 g。盡量選擇在平穩(wěn)耕地、農(nóng)用菜地、林果地、山坡下側(cè)土層較厚等地帶采樣，避開明顯點(diǎn)狀污染或表土已被破壞的干擾地段，保證了樣品的典型性和代表性。

測(cè)試的指標(biāo)為類重金屬As，重金屬元素Cd、Cr、Cu、Hg、Ni、Pb、Zn，微量元素Se。采用氫化物發(fā)生原子熒光儀測(cè)定As和Se，采用ICP-OES (PE，USA)測(cè)定Cd、Cr、Cu、Hg、Ni、Pb、Zn。分析測(cè)試中準(zhǔn)確度和精密度采用國(guó)家一級(jí)土壤標(biāo)準(zhǔn)物質(zhì)(GBW07349)控制，加10%空白樣與平行樣控制，且合格率符合規(guī)范要求，指標(biāo)的加標(biāo)回收率均符合國(guó)家標(biāo)準(zhǔn)。

2.2 數(shù)據(jù)處理

應(yīng)用Excel 2016進(jìn)行土壤數(shù)據(jù)基本參數(shù)統(tǒng)計(jì)分析，應(yīng)用SPSS 19.0完成土壤數(shù)據(jù)正態(tài)分布性檢驗(yàn)、相關(guān)系數(shù)分析、聚類分析、因子分析等。為了避免在利用統(tǒng)計(jì)學(xué)和地統(tǒng)計(jì)學(xué)方法時(shí)可能出現(xiàn)的比例效應(yīng)[14]，所分析的數(shù)據(jù)必須全部服從正態(tài)分布或近似正態(tài)分布，否則需要將數(shù)據(jù)進(jìn)行對(duì)數(shù)變換或者Box-Cox變換。聚類分析方法采用最近鄰元素，選擇區(qū)間Pearson相關(guān)性作為度量標(biāo)準(zhǔn)。因子分析采用最大四次方值法，抽取基于特征值大于1，最大收斂性迭代次數(shù)為25次。

目前國(guó)內(nèi)外判斷數(shù)據(jù)是否服從正態(tài)分布的方法較多[15]，主要與所使用的統(tǒng)計(jì)學(xué)軟件類型(SAS、STATA、SPSS等)和樣本量的多少有關(guān)。使用SPSS軟件時(shí)，規(guī)則為探索性檢驗(yàn)中基于理論假定的概率圖、Q-Q圖[16]、P-P圖等不限定樣本量；偏度峰度檢驗(yàn)法適用于樣本量N>200時(shí)；樣本量72 000時(shí)，以Lilliefors檢驗(yàn)(K-S檢驗(yàn)[17])為準(zhǔn)?；谔幚碥浖蜆颖緮?shù)量的限定，本文選用偏度峰度檢驗(yàn)法和W檢驗(yàn)法，結(jié)果表明，偏度值和峰度值遠(yuǎn)大于1，W檢驗(yàn)的顯著性均小于0.01；對(duì)9項(xiàng)元素進(jìn)行對(duì)數(shù)變換后偏度值和峰度值基本在1附近，W檢驗(yàn)的顯著性均大于0.05，符合正態(tài)分布的要求，滿足地統(tǒng)計(jì)學(xué)分析的假設(shè)條件，可以進(jìn)行半變異函數(shù)的計(jì)算。

應(yīng)用ArcGIS 10.2中Geostatistics analysis模塊的普通克里格法結(jié)合半方差模型，完成最優(yōu)內(nèi)插過程。

應(yīng)用GS+9.0軟件完成9項(xiàng)元素半變異函數(shù)的計(jì)算和多種理論模型的擬合。半變異函數(shù)作為描述隨機(jī)場(chǎng)和隨機(jī)過程空間相關(guān)性的重要統(tǒng)計(jì)量，是分析土壤指標(biāo)空間變異結(jié)構(gòu)性和隨機(jī)性的有效性方法[18]。該函數(shù)如下：

(1)

式中：γ(h)——半變異函數(shù)值；

N(h)——樣點(diǎn)對(duì)的數(shù)目；

h——步長(zhǎng)，即兩分隔樣點(diǎn)的距離；

Z(xi)——Z(xi)在空間位置xi上的土壤觀測(cè)值；

Z(xi+h)——Z(xi)在空間位置xi+h上的土壤觀測(cè)值。

在土壤地球化學(xué)工作中，通常區(qū)域化變量往往在不同的方向上會(huì)表現(xiàn)出不同的變異性，即便是在同一方向上也存在不同尺度的多維度的變異性。為了精準(zhǔn)反映區(qū)域化變量的主要空間變異性，需要根據(jù)半變異函數(shù)的決定系數(shù)(R2)和殘差(RSS)按照土壤元素所對(duì)應(yīng)的最優(yōu)理論模型進(jìn)行半變異函數(shù)的擬合。本次所涉及的半方差函數(shù)模型主要有以下兩種：

指數(shù)模型(Exponential)：

(2)

線性模型(Linear)：

γ(h)=C0+C·h/α，h>0

(3)

式中：C0——塊金常數(shù)(間距為0時(shí)的半方差，通常由隨機(jī)因素引起的變異)；

C——拱高(通常由系統(tǒng)因素引起的變異)；

α——變程(半方差達(dá)到基臺(tái)值的樣本間距，體現(xiàn)隨機(jī)變量在空間上的自相關(guān)程度)。

3 結(jié)果分析

3.1 土壤重金屬元素的統(tǒng)計(jì)分析

研究區(qū)內(nèi)共采集土壤表層樣品401件，表層土壤重金屬元素是以統(tǒng)計(jì)學(xué)為基本依據(jù)，通過系統(tǒng)性計(jì)算各指標(biāo)的最小值、最大值、中位數(shù)、平均值、標(biāo)準(zhǔn)差等統(tǒng)計(jì)特征值，來描述表層土壤重金屬元素的含量區(qū)間和分布規(guī)律，見表1。表層土壤重金屬元素的中位數(shù)與平均值相差較大，表明其中心趨向分布可能被異常值影響而使其呈非標(biāo)準(zhǔn)正態(tài)分布。將各指標(biāo)平均值與河北省背景值[19]對(duì)比，發(fā)現(xiàn)土壤中Cd、Cu和Pb的平均含量均高于河北省背景值，其中以Cd最高，是河北省背景值的2倍；As、Cr、Hg、Ni和Zn的平均含量均與河北省背景值基本一致。

表1 土壤重金屬的描述性統(tǒng)計(jì)

大量研究表明，變異系數(shù)(Cv)不僅可以反映總體土壤樣本中各樣點(diǎn)元素含量的平均變異程度，也可反映土壤元素的離散程度[20]。變異系數(shù)值越大，元素的空間分布越不均勻，離散程度越高。參照本次測(cè)試數(shù)據(jù)特點(diǎn)，將變異程度劃分為三種類型：Cv<70%為均勻分布、70%≤Cv<100%為中強(qiáng)分異(中高起伏)、Cv≥100%為強(qiáng)分異(很大起伏)。依據(jù)該分類標(biāo)準(zhǔn)，Cd、Cu、Hg和Pb的變異系數(shù)分別為385%、143%、350%、118%，表明表層土壤中Cd、Cu、Hg和Pb屬于強(qiáng)分異，其中以Cd和Hg的分異程度最強(qiáng)，表明局部地區(qū)可能受到人為活動(dòng)影響，存在點(diǎn)源污染；As、Cr、Ni的變異系數(shù)分別為96%、93%、86%，屬于中強(qiáng)分異；Zn的變異系數(shù)為64%，分異程度最弱。

在置信區(qū)間為95%內(nèi)，不同土壤類型中單因素方差分析結(jié)果為F=2.913、p=0.006，判斷土壤類型對(duì)土壤Zn含量均值具有顯著性影響(表2)；不同土地利用類型中單因素方差分析結(jié)果表明，Cr的F=3.408、p=0.003，Cu的F=3.097、p=0.006，Ni的F=2.355、p=0.030，判斷土地利用類型對(duì)土壤Cr、Cu、Ni含量均值具有顯著性影響(表3)。

因此按照不同土壤類型，對(duì)Zn含量均值進(jìn)行統(tǒng)計(jì)分析，發(fā)現(xiàn)潮土>栗鈣土>粗骨土>褐土>棕壤>草甸土>灰色森林土>山地草甸土。按照不同土地利用類型，對(duì)Cr、Cu、Ni含量均值進(jìn)行統(tǒng)計(jì)分析，發(fā)現(xiàn)Cu在內(nèi)陸灘涂中含量高，Cr和Ni在建設(shè)用地中含量高，可能與人為活動(dòng)影響有關(guān)。

表2 9種元素在不同土壤類型中的均值特征及單方差分析結(jié)果

表3 9種元素在不同土地利用類型的均值特征及單方差分析結(jié)果

3.2 土壤元素之間的相關(guān)性分析

對(duì)土壤重金屬元素和Se進(jìn)行相關(guān)性分析(spearman分析)，相關(guān)系數(shù)R越大，元素之間的關(guān)聯(lián)程度越高；相關(guān)系數(shù)R越小，元素之間的關(guān)聯(lián)程度越低或同源性越差。由各元素間的相關(guān)系數(shù)表可知(表4)，具有顯著性相關(guān)關(guān)系的約占總關(guān)系數(shù)的25%，極顯著相關(guān)關(guān)系約占總關(guān)系數(shù)的47.2%，其中Cr和Ni、Cd和Pb、Cd和Zn、Pb和Zn、Cu和Hg呈極顯著的正相關(guān)，相關(guān)系數(shù)分別為0.832，0.698，0.739，0.731，0.635。Se與As、Hg、Pb、Zn呈極顯著的正相關(guān)，與Cu和Ni呈顯著的正相關(guān)，而與Cd和Cr之間的相關(guān)性較弱。Cr與As和Pb之間呈一定的負(fù)相關(guān)關(guān)系。由此推斷Cr和Ni、Cd和Pb和Zn、Cu和Hg、Se和As之間可能具有相同的來源。

表4 9個(gè)土壤元素的相關(guān)性分析

3.3 土壤元素聚類及主成分分析

通過對(duì)研究區(qū)表層土壤中9種元素?cái)?shù)據(jù)進(jìn)行聚類分析，樹狀圖結(jié)果表明，Cr和Ni最先聚合為一類，隨后Cd、Pb、Zn聚合為第一類，Cu和Hg聚合為第二類，而As和Se聚合為第三類，如圖2所示。

在一定程度上通常采用特征值來表征主成分影響力度大小，本次提取特征值大于1的主成分。由主成分分析結(jié)果可知(表5)，表層土壤中9種元素可被劃分為四個(gè)主成分，第一個(gè)主成分(F1)的方差貢獻(xiàn)率為34.50%，占四種主成分中較大比例，可以作為衡量該區(qū)重金屬狀況的一個(gè)綜合性指標(biāo)。第二個(gè)主成分(F2)的方差貢獻(xiàn)率為20.55%，第三個(gè)主成分(F3)的方差貢獻(xiàn)率為14.25%，第四個(gè)主成分(F4)的方差貢獻(xiàn)率為12.58%，這四個(gè)成分的累積方差貢獻(xiàn)率為81.87%，具備原始數(shù)據(jù)的大部分信息。

表5 土壤元素的主成分分析結(jié)果

為了更加精準(zhǔn)確定各主成分所包含的元素信息，經(jīng)過方差最大正交旋轉(zhuǎn)后，對(duì)9種元素的四種主成分上因子載荷進(jìn)行統(tǒng)計(jì)，如圖3。結(jié)果發(fā)現(xiàn)，Cd、Pb、Zn在第一主成分中為高載荷，而在第二主成分、第三主成分、第四主成分中均為低載荷；Cr和Ni在第二主成分中為高載荷，而在第一主成分、第三主成分、第四主成分中均為低載荷；Cu和Hg在第三主成分中為高載荷，而在第一主成分、第二主成分、第四主成分中均為低載荷；As和Se在第四主成分中為高載荷，而在第一主成分、第二主成分、第三主成分中均為低載荷，見表6。

圖3 表層土壤元素主成分載荷Fig.3 Factor loading analysis of elements in surficial soils

綜上所述，傳統(tǒng)統(tǒng)計(jì)學(xué)中主成分分析與相關(guān)性分析和聚類分析的結(jié)果在一定程度上保持一致，即Cr和Ni，Cd、Pb、Zn三者，Cu和Hg，As和Se來源可能相同。

表6 主成分分析旋轉(zhuǎn)后的成分載荷矩陣

3.4 土壤元素的趨勢(shì)分析

表層土壤一般容易受到外界因素(氣候條件、地形地貌、人為活動(dòng)等[21])的共同影響，因此區(qū)域性的表層土壤空間分布具有顯著的趨勢(shì)性和異向性。目前可通過Arcgis 10.2中趨勢(shì)分析生成的三維趨勢(shì)面曲線形狀以及其變化的陡峭程度來判斷空間數(shù)據(jù)的總體變化趨勢(shì)。一般分為三種類型，即無(沒有趨勢(shì)效應(yīng))、一階(區(qū)域化變量沿一定方向呈直線變化)和二階(區(qū)域化變量沿一定方向呈多項(xiàng)式變化)[22]，通過異向性軸向自動(dòng)搜索模塊，對(duì)各元素?cái)?shù)據(jù)進(jìn)行處理，從而確定各向異性的趨勢(shì)效應(yīng)特征。以研究區(qū)As、Cu、Cd的空間分布趨勢(shì)效應(yīng)分析圖為例，分別代表了無、一階、二階趨勢(shì)效應(yīng)(圖4)，圖中東西方向由x軸表示，南北方向由y軸表示，每個(gè)樣點(diǎn)實(shí)測(cè)值的大小由z軸表示；正北方向投影面上淺綠色曲線表示東西向的全局性趨勢(shì)效應(yīng)變化，正東向投影面上紫色曲線表示南北向的全局性趨勢(shì)效應(yīng)變化。從圖上可以看出，研究區(qū)土壤中Cd表現(xiàn)為在東西方向上先減小后增大的趨勢(shì)，Cu表現(xiàn)為在東西方向和南北方向上呈直線分布，而As表現(xiàn)為無階效應(yīng)。

圖4 土壤元素趨勢(shì)分析圖Fig.4 Trend analysis of soil elements

在Kringing插值過程中，選用不同的趨勢(shì)類型會(huì)導(dǎo)致插值過程中產(chǎn)生的誤差不同，因此需要綜合分析各類參數(shù)特點(diǎn)，最大程度將各類誤差降為最低，進(jìn)而使半方差函數(shù)模型及其參數(shù)最適合。本次研究參照以下標(biāo)準(zhǔn)判斷，平均誤差(mean error，ME)的絕對(duì)值最接近0，均方根誤差(root-mean-square error，RMSE)與平均標(biāo)準(zhǔn)誤差(average standard error，ASE)應(yīng)盡量相等且值越小越優(yōu)，若前者小于后者，表明低估了預(yù)測(cè)值，反之則說明高估了預(yù)測(cè)值。標(biāo)準(zhǔn)化均方根誤差(root-mean-square standardized error，RMSSE)值最接近1，如果標(biāo)準(zhǔn)化均方根誤差小于1，說明高估了預(yù)測(cè)值，反之為低估預(yù)測(cè)值[23]，標(biāo)準(zhǔn)化平均誤差(mean standardized error，MSE)的絕對(duì)值應(yīng)最接近0。由表7可知，As、Cr、Pb、Zn、Ni、Se的無趨勢(shì)預(yù)測(cè)與一階趨勢(shì)預(yù)測(cè)和二階趨勢(shì)預(yù)測(cè)相比，ME的絕對(duì)值更加接近0，ASE和RMSE最接近，MSE的絕對(duì)值更加接近0，RMSSE的值最接近1，因此Kringing插值時(shí)應(yīng)選擇無趨勢(shì)更加精準(zhǔn)，同理Hg和Cu選擇一階，Cd選擇二階。

3.5 土壤元素空間異質(zhì)性特征

半變異函數(shù)(Semivariogram)是在區(qū)域化變量滿足二階平穩(wěn)和本征假設(shè)的前提下，按照一定的抽樣間隔獲得樣本方差的數(shù)學(xué)期望值，是地統(tǒng)計(jì)學(xué)中分析空間格局變量時(shí)最主要的應(yīng)用工具。通常應(yīng)用半變異函數(shù)及其參數(shù)可以對(duì)區(qū)域化變量的分布進(jìn)行結(jié)構(gòu)性和隨機(jī)性的判斷，進(jìn)而結(jié)合克里格法完成對(duì)無樣品空白區(qū)的插值和預(yù)測(cè)[24]。有三個(gè)相應(yīng)的參數(shù)可直接反應(yīng)空間變異性程度，塊金值(C0)是在取樣距離為0時(shí)，半變異函數(shù)為一個(gè)非0定值，反映的是最小抽樣尺度以下由儀器測(cè)量誤差等隨機(jī)因素引起的微變異；變程(A0)也稱半變異函數(shù)達(dá)到基臺(tái)值時(shí)空間最大間隔距離，反映了變量空間自相關(guān)范圍；塊金系數(shù)C0/(C+C0)是塊金值與基臺(tái)值的比值，反應(yīng)了空間異質(zhì)性。而空間變異性程度主要由隨機(jī)性變異和結(jié)構(gòu)性變異大小所決定，結(jié)構(gòu)性因素會(huì)增加變量的空間相關(guān)性，隨機(jī)性因素則會(huì)降低變量的空間相關(guān)性，當(dāng)塊金系數(shù)小于25%時(shí)，表明區(qū)域變量的空間相關(guān)性強(qiáng)烈，此時(shí)空間變異主要由地形地貌、土壤類型、土地利用方式等結(jié)構(gòu)性因素控制[25～26]；若塊金系數(shù)大于75%，則表明區(qū)域變量幾乎不具有或具有較弱的空間相關(guān)性，其空間變異程度受隨機(jī)性因子控制比重較高；當(dāng)塊金系數(shù)介于25%～75%，表明區(qū)域變量的空間相關(guān)性強(qiáng)度居中，其空間變異由隨機(jī)性因素與結(jié)構(gòu)性因素共同控制[27]。土壤表層重金屬元素空間變異性研究的關(guān)鍵是半變異函數(shù)模型的確定，模型擬合的效果可按照R2值最大和RSS值最小的原則檢驗(yàn)[28]。

由表8可知，9種土壤元素的R2值除了Pb為0.176，其余均在0.334～0.773之間，說明上述重金屬元素和Se的理論半變異與實(shí)驗(yàn)半變異模型較為吻合，擬合效果趨優(yōu)。這9種元素中有8項(xiàng)元素的變程遠(yuǎn)遠(yuǎn)超過了8 000 m，最小的變程也在3 000 m以上，因此可對(duì)其進(jìn)行Kriging插值。

As的半方差函數(shù)屬于指數(shù)模型，塊金系數(shù)為10.8%，變程為3 000 m，因此As具強(qiáng)烈的空間相關(guān)性，反映出研究區(qū)表層土壤中As主要受到地形地貌、成土母質(zhì)、地質(zhì)背景等結(jié)構(gòu)性因素的影響，這與Chen等[29]和何廳廳等[30]結(jié)論保持一致，他們認(rèn)為As的空間分布主控因素為土壤屬性和地質(zhì)背景，受人為活動(dòng)影響較??；Cr、Cd、Hg、Pb、Zn、Cu、Ni、Se的塊金系數(shù)分別為83.7%、83.0%、87.4%、100%、81.7%、89.6%、85.7%、86.9%，且半方差函數(shù)均屬于線性模型(圖5)。

表7 不同趨勢(shì)階數(shù)插值誤差比較

表8 土壤元素含量變異函數(shù)理論模型及其相關(guān)參數(shù)

這8項(xiàng)元素具較弱的空間相關(guān)性，表明研究區(qū)表層土壤中這8項(xiàng)元素主要受到人為活動(dòng)等隨機(jī)性因素的影響，結(jié)論與Facchinelli等的[31]一致，其中Cu、Zn和Pb的來源與人類活動(dòng)較為密切。不存在中等強(qiáng)度的空間變異性情況，這可能也與本次研究區(qū)邊緣地帶樣品點(diǎn)較為稀疏有一定關(guān)系。

圖5 土壤元素半方差函數(shù)圖Fig.5 Semivariogram of the soil elements

3.6 普通克里格(Kriging)空間插值法

GIS在處理海量數(shù)據(jù)和圖形視覺化表達(dá)方面具有很突出的特點(diǎn)，選取模型插值成圖后可以直觀地反映土壤元素的空間變異性?？臻g插值法能夠使用一定量的樣本量對(duì)區(qū)域范圍內(nèi)空白樣本進(jìn)行預(yù)測(cè)，從而分析整個(gè)區(qū)域范圍的相關(guān)特征?？臻g插值的方法較多，主要分為確定性插值法(反距離加權(quán)插值法、徑向基函數(shù)插值法、局部多項(xiàng)式插值等[32-33])和地統(tǒng)計(jì)學(xué)插值法(克里格法/協(xié)同克里格法、面插值法、經(jīng)驗(yàn)貝葉斯克里金法等[34-35])?？死锔穹ㄓ址Q空間局部插值法，是在一定區(qū)域內(nèi)對(duì)區(qū)域化變量進(jìn)行無偏最優(yōu)估計(jì)的一種方法，可分為泛克里格插值、簡(jiǎn)單克里格插值等，由于滿足內(nèi)蘊(yùn)假設(shè)等條件[36]，本文采用地統(tǒng)計(jì)學(xué)中常用的普通克里格空間插值法[37]。

圖6 土壤元素的空間分布Fig.6 Spatial distribution of soil elements

選取最優(yōu)趨勢(shì)參數(shù)和半方差函數(shù)模型進(jìn)行普通克里格空間插值，得到研究區(qū)土壤元素的空間分布圖(圖6)。參照《土地質(zhì)量地球化學(xué)評(píng)價(jià)規(guī)范》[38]和河北省土壤元素背景值調(diào)查[19]， As含量范圍為2.9～29.9 mg/kg，整體含量值偏低；Ni含量范圍為12.7～61.9 mg/kg，Pb含量范圍為16.4～38.6 mg/kg，Se含量范圍為0.055～0.401 mg/kg，整體含量值處于中等；Cd含量范圍為0.082～0.742 mg/kg，Cr含量范圍為28.0～228.4 mg/kg，Cu含量范圍為11.3～142.4 mg/kg，Hg含量范圍0.013～0.126 mg/kg，Zn含量范圍為43.5～199.4 mg/kg，整體含量值偏高。從區(qū)域空間分布來看，承德市北部地區(qū)9種元素含量整體偏低，而南部地區(qū)含量整體偏高，體現(xiàn)出壩上高原地區(qū)整體受人為活動(dòng)影響較低的特點(diǎn)。Cr和Ni空間分布具有很強(qiáng)的一致性，以豐寧縣和隆化縣為界，低值區(qū)主要分布于承德地區(qū)的北部，南部鷹手營(yíng)子礦區(qū)和寬城滿族自治縣附近有低值區(qū)的零星分布；Cu和Hg保持了一定的空間分布一致性，高值區(qū)主要分布于承德南部地區(qū)，而Hg在豐寧縣北部地區(qū)亦有高值區(qū)；Pb在中部形成了一條高值帶，這又與Cd分布特征相似，這與區(qū)內(nèi)鐵礦、金礦等金屬礦床的集中分布狀況具有一定相關(guān)性。除去Pb中部的高值帶，Zn和Pb二者在承德南部以及中部地區(qū)的高值區(qū)近似分布， Ni的高值區(qū)較好地表現(xiàn)出由南向北逐漸遞變的過度特性，可能與區(qū)內(nèi)金屬礦產(chǎn)分布狀況有一定關(guān)系；As低值區(qū)主要分布于承德中部隆化縣—灤平縣地區(qū)，圍場(chǎng)滿族蒙古族自治縣以北、興隆縣以西地區(qū)。雖然上文統(tǒng)計(jì)學(xué)分析并未提及Se和Cd具有某種同源關(guān)系，但是二者高值區(qū)都分布在承德南部地區(qū)，西北部地區(qū)亦呈現(xiàn)局部高值，灤平縣—雙灤區(qū)、寬城滿族自治縣達(dá)到最高值，而西南部邊緣地帶和北部地區(qū)呈現(xiàn)低值區(qū)，也表現(xiàn)出弱空間異構(gòu)性，這可能與Cd和Se都具有親硫性和親生物性，土壤中Cd、Se經(jīng)常共生在一起有關(guān)，因此為富硒土地資源的安全開發(fā)利用方面帶來了嚴(yán)峻挑戰(zhàn)。

綜上所述，單憑某一種方法去判斷元素之間的同源關(guān)系是不準(zhǔn)確的，需要在傳統(tǒng)統(tǒng)計(jì)學(xué)結(jié)論的基礎(chǔ)上進(jìn)一步參考變異函數(shù)模型并疊加克里格插值結(jié)果進(jìn)行優(yōu)化。綜合分析認(rèn)為，研究區(qū)Cr和Ni、Cu和Hg、Zn和Pb、Cd和Se之間具有很好的空間分布相關(guān)性，來源具有一致性。建議應(yīng)在現(xiàn)有土地利用現(xiàn)狀屬性基礎(chǔ)上，充分考慮土壤重金屬元素空間分布格局要素，在有礦業(yè)活動(dòng)且土壤呈現(xiàn)重金屬高值地區(qū)(如灤平縣—雙灤區(qū)北部山區(qū))進(jìn)行周期性土壤監(jiān)測(cè)。同時(shí)，進(jìn)一步開展區(qū)內(nèi)土壤Cd和Se兩者定量化控制因素的厘定及模型機(jī)理的研究，從而實(shí)現(xiàn)農(nóng)作物富硒降鎘的目標(biāo)。

4 結(jié)論

(1)傳統(tǒng)統(tǒng)計(jì)學(xué)分析中Cd、Cu、Hg和Pb屬于強(qiáng)分異，變異系數(shù)分別為385%、143%、350%、118%，As、Cr和Ni屬于中強(qiáng)分異，變異系數(shù)分別為96%、93%、86%，Zn屬于弱分異，變異系數(shù)為64%。Zn含量均值受土壤類型影響顯著，Cr、Cu、Ni含量均值則受土地利用類型影響顯著。結(jié)合相關(guān)性分析、聚類分析和主成分結(jié)果，可知Cr和Ni、Cd和Pb和Zn、Cu和Hg來源相同。

(2)通過異向性軸向自動(dòng)搜索功能結(jié)合不同趨勢(shì)階數(shù)插值誤差綜合對(duì)比，確定了As、Cr、Pb、Zn、Ni、Se適宜選擇無趨勢(shì)，Hg和Cu適宜選擇一階，而Cd適宜選擇二階。

(3)地統(tǒng)計(jì)學(xué)分析結(jié)果表明，As采用指數(shù)模型可較好地?cái)M合，Cd、Cr、Cu、Hg、Ni、Pb、Zn、Se采用線性模型可較好地?cái)M合。其中As的塊金效應(yīng)小于25%，因此具強(qiáng)烈的空間相關(guān)性，受結(jié)構(gòu)性因素影響較大；Cd、Cr、Cu、Hg、Ni、Pb、Zn、Se的塊金效應(yīng)大于75%，具較弱的空間相關(guān)性，受隨機(jī)性因素影響較大。

(4)普通克里金插值圖直觀地反映了承德地區(qū)土壤重金屬的空間結(jié)構(gòu)特征。研究區(qū)土壤的9種重金屬元素都呈現(xiàn)明顯的北低南高趨勢(shì)，Cr和Ni、Cu和Hg、Zn和Pb、Se和Cd之間高值區(qū)分布具有很強(qiáng)的空間一致性，Hg在豐寧縣北部地區(qū)亦有高值區(qū)，Zn和Pb高值區(qū)都分布于承德南部有向中部地區(qū)發(fā)展的趨勢(shì)， Pb在中部地區(qū)形成了一條較寬的高值帶，而其余分布特征與Cd相似，Ni的高值區(qū)較好地表現(xiàn)出由南向北逐漸遞變的過度特性，在一定程度上較好地印證了地統(tǒng)計(jì)學(xué)所得出的規(guī)律。