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

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

土壤是成土母巖經(jīng)過風(fēng)化、剝蝕、搬運(yùn)、沉積等地質(zhì)作用后疊加成土作用形成的自然物質(zhì)體。在成土的每個(gè)階段，都可能會受到外界環(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]。目前已有不少國內(nèi)外學(xué)者針對不同區(qū)域時(shí)空尺度，利用地統(tǒng)計(jì)學(xué)結(jié)合GIS的方法解決土壤元素空間變異性問題。如崔萌等[5]運(yùn)用地統(tǒng)計(jì)學(xué)和主成分分析法對北京市大桃主產(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é)方法對三峽庫區(qū)表層土壤部分微量元素進(jìn)行了空間變異特征分析，結(jié)果表明7種微量元素為中等自相關(guān)性，2種為弱自相關(guān)性，且4個(gè)方向上變異程度不明顯；Roger等[7]通過環(huán)境因子作為數(shù)據(jù)對土壤養(yǎng)分指數(shù)完成克里格插值，表明克里格插值與線性回歸共同分析可以一定程度提高土壤養(yǎng)分的預(yù)測精度。

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

1 研究區(qū)概況

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

承德礦產(chǎn)資源分布較多，以鐵礦資源最為豐富，是我國僅次于攀枝花的第二大釩鈦磁鐵礦資源基地[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)狀污染或表土已被破壞的干擾地段，保證了樣品的典型性和代表性。

測試的指標(biāo)為類重金屬As，重金屬元素Cd、Cr、Cu、Hg、Ni、Pb、Zn，微量元素Se。采用氫化物發(fā)生原子熒光儀測定As和Se，采用ICP-OES (PE，USA)測定Cd、Cr、Cu、Hg、Ni、Pb、Zn。分析測試中準(zhǔn)確度和精密度采用國家一級土壤標(biāo)準(zhǔn)物質(zhì)(GBW07349)控制，加10%空白樣與平行樣控制，且合格率符合規(guī)范要求，指標(biāo)的加標(biāo)回收率均符合國家標(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)行對數(shù)變換或者Box-Cox變換。聚類分析方法采用最近鄰元素，選擇區(qū)間Pearson相關(guān)性作為度量標(biāo)準(zhǔn)。因子分析采用最大四次方值法，抽取基于特征值大于1，最大收斂性迭代次數(shù)為25次。

目前國內(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；對9項(xiàng)元素進(jìn)行對數(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ī)場和隨機(jī)過程空間相關(guān)性的重要統(tǒng)計(jì)量，是分析土壤指標(biāo)空間變異結(jié)構(gòu)性和隨機(jī)性的有效性方法[18]。該函數(shù)如下：

(1)

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

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

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

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

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

在土壤地球化學(xué)工作中，通常區(qū)域化變量往往在不同的方向上會表現(xiàn)出不同的變異性，即便是在同一方向上也存在不同尺度的多維度的變異性。為了精準(zhǔn)反映區(qū)域化變量的主要空間變異性，需要根據(jù)半變異函數(shù)的決定系數(shù)(R2)和殘差(RSS)按照土壤元素所對應(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á)到基臺值的樣本間距，體現(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]對比，發(fā)現(xiàn)土壤中Cd、Cu和Pb的平均含量均高于河北省背景值，其中以Cd最高，是河北省背景值的2倍；As、Cr、Hg、Ni和Zn的平均含量均與河北省背景值基本一致。

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

大量研究表明，變異系數(shù)(Cv)不僅可以反映總體土壤樣本中各樣點(diǎn)元素含量的平均變異程度，也可反映土壤元素的離散程度[20]。變異系數(shù)值越大，元素的空間分布越不均勻，離散程度越高。參照本次測試數(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，判斷土壤類型對土壤Zn含量均值具有顯著性影響(表2)；不同土地利用類型中單因素方差分析結(jié)果表明，Cr的F=3.408、p=0.003，Cu的F=3.097、p=0.006，Ni的F=2.355、p=0.030，判斷土地利用類型對土壤Cr、Cu、Ni含量均值具有顯著性影響(表3)。

因此按照不同土壤類型，對Zn含量均值進(jìn)行統(tǒng)計(jì)分析，發(fā)現(xiàn)潮土>栗鈣土>粗骨土>褐土>棕壤>草甸土>灰色森林土>山地草甸土。按照不同土地利用類型，對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)性分析

對土壤重金屬元素和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 土壤元素聚類及主成分分析

通過對研究區(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)后，對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 土壤元素的趨勢分析

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

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

在Kringing插值過程中，選用不同的趨勢類型會導(dǎo)致插值過程中產(chǎn)生的誤差不同，因此需要綜合分析各類參數(shù)特點(diǎn)，最大程度將各類誤差降為最低，進(jìn)而使半方差函數(shù)模型及其參數(shù)最適合。本次研究參照以下標(biāo)準(zhǔn)判斷，平均誤差(mean error，ME)的絕對值最接近0，均方根誤差(root-mean-square error，RMSE)與平均標(biāo)準(zhǔn)誤差(average standard error，ASE)應(yīng)盡量相等且值越小越優(yōu)，若前者小于后者，表明低估了預(yù)測值，反之則說明高估了預(yù)測值。標(biāo)準(zhǔn)化均方根誤差(root-mean-square standardized error，RMSSE)值最接近1，如果標(biāo)準(zhǔn)化均方根誤差小于1，說明高估了預(yù)測值，反之為低估預(yù)測值[23]，標(biāo)準(zhǔn)化平均誤差(mean standardized error，MSE)的絕對值應(yīng)最接近0。由表7可知，As、Cr、Pb、Zn、Ni、Se的無趨勢預(yù)測與一階趨勢預(yù)測和二階趨勢預(yù)測相比，ME的絕對值更加接近0，ASE和RMSE最接近，MSE的絕對值更加接近0，RMSSE的值最接近1，因此Kringing插值時(shí)應(yīng)選擇無趨勢更加精準(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ù)可以對區(qū)域化變量的分布進(jìn)行結(jié)構(gòu)性和隨機(jī)性的判斷，進(jìn)而結(jié)合克里格法完成對無樣品空白區(qū)的插值和預(yù)測[24]。有三個(gè)相應(yīng)的參數(shù)可直接反應(yīng)空間變異性程度，塊金值(C0)是在取樣距離為0時(shí)，半變異函數(shù)為一個(gè)非0定值，反映的是最小抽樣尺度以下由儀器測量誤差等隨機(jī)因素引起的微變異；變程(A0)也稱半變異函數(shù)達(dá)到基臺值時(shí)空間最大間隔距離，反映了變量空間自相關(guān)范圍；塊金系數(shù)C0/(C+C0)是塊金值與基臺值的比值，反應(yīng)了空間異質(zhì)性。而空間變異性程度主要由隨機(jī)性變異和結(jié)構(gòu)性變異大小所決定，結(jié)構(gòu)性因素會增加變量的空間相關(guān)性，隨機(jī)性因素則會降低變量的空間相關(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以上，因此可對其進(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)影響較?。籆r、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ù)插值誤差比較

表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插值法能夠使用一定量的樣本量對區(qū)域范圍內(nèi)空白樣本進(jìn)行預(yù)測，從而分析整個(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)對區(qū)域化變量進(jìn)行無偏最優(yōu)估計(jì)的一種方法，可分為泛克里格插值、簡單克里格插值等，由于滿足內(nèi)蘊(yùn)假設(shè)等條件[36]，本文采用地統(tǒng)計(jì)學(xué)中常用的普通克里格空間插值法[37]。

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

選取最優(yōu)趨勢參數(shù)和半方差函數(shù)模型進(jìn)行普通克里格空間插值，得到研究區(qū)土壤元素的空間分布圖(圖6)。參照《土地質(zhì)量地球化學(xué)評價(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ū)的北部，南部鷹手營子礦區(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ū)，圍場滿族蒙古族自治縣以北、興隆縣以西地區(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)測。同時(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ù)插值誤差綜合對比，確定了As、Cr、Pb、Zn、Ni、Se適宜選擇無趨勢，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)明顯的北低南高趨勢，Cr和Ni、Cu和Hg、Zn和Pb、Se和Cd之間高值區(qū)分布具有很強(qiáng)的空間一致性，Hg在豐寧縣北部地區(qū)亦有高值區(qū)，Zn和Pb高值區(qū)都分布于承德南部有向中部地區(qū)發(fā)展的趨勢， Pb在中部地區(qū)形成了一條較寬的高值帶，而其余分布特征與Cd相似，Ni的高值區(qū)較好地表現(xiàn)出由南向北逐漸遞變的過度特性，在一定程度上較好地印證了地統(tǒng)計(jì)學(xué)所得出的規(guī)律。