徐克全 ， 金立新， 鄭勇軍， 張 華， 包雨函 ， 時章亮

(1.四川省地質(zhì)調(diào)查院，成都 610081；2.稀有稀土戰(zhàn)略資源評價與利用四川省重點實驗室，成都 610081)

0 引言

土壤元素含量在空間上具有自相關(guān)性，也存在多元素間交叉相關(guān)[1-2]，在空間分布上總體受到地質(zhì)背景的控制，局部可能會受到人為因素的影響。為了獲取特定空間范圍數(shù)據(jù)，并對特定空間對象進行評價，需要對有限的空間數(shù)據(jù)進行插值處理。常用到反距離加權(quán)(IDW)、半變異函數(shù)的克里金(Kriging)等方法進行插值。但不同方法選取，或半變異函數(shù)中的塊金、變程、偏基臺值的設置，都將對插值結(jié)果造成影響，最終影響土壤質(zhì)量地球化學綜合評價面積統(tǒng)計的準確性，不利于土地風險管控和科學利用。筆者通過R語言(運用到sp、gstat、geoR、ggplot2、raster、intensity.analysis等R語言包)，使用麻城鎮(zhèn)和摩尼鎮(zhèn)1:10 000土地質(zhì)量地球化學調(diào)查評價結(jié)果的表層土壤數(shù)據(jù)[3]，考慮元素空間自相關(guān)性可能受到地質(zhì)背景、土地利用類型、元素的地球化學親和性、空間距離的影響，評價單元對應賦值方法分別采用地層屬性分層克里金法插值、土地利用類型分層克里金法插值、多變量協(xié)同克里金法插值、反距離加權(quán)平均插值法；考慮“有樣”評價單元賦原值以及“無樣”評價單元插值方法的更合理性，也是主流方法，增加對“有樣”單元均采用局部加權(quán)平均(權(quán)重為公式3)賦值，而對“無樣”單元插值分別采用“地層屬性分層克里金插值”、“土地利用類型分層克里金插值”、“多變量協(xié)同克里金插值”和“反距離加權(quán)平均法(權(quán)重公式(2))插值”，共8種賦值方法進行養(yǎng)分單指標和土壤質(zhì)量地球化學綜合評價。借鑒土地利用強度分析方法，對以上8種賦值結(jié)果進行對比分析，探討各種方法對評價結(jié)果的差異。

1 研究區(qū)概況

研究區(qū)位于四川省敘永縣東南部麻城鎮(zhèn)與摩尼鎮(zhèn)境內(nèi)，面積為54.59 km2，其中麻城鎮(zhèn)于2015年撤鄉(xiāng)設鎮(zhèn)(圖1)。共采集1 791(含重復樣40件)件樣品，基本密度為32.07件/km2；重復樣數(shù)為40件，占2.29%；按網(wǎng)格化方法進行采樣，采樣單元以0.25 km2為采樣大格，250 m×250 m采樣小格，以采樣小格作為評價單元。

圖1 研究區(qū)地理位置圖[3]

出露地層復雜多樣，主要有寒武系(ε)、奧陶系(O)、志留系(S)、二疊系(P)和第四系(Q)(圖2)。

圖2 研究區(qū)地質(zhì)圖[3]

1)寒武系。婁山關(guān)組(∈2-O1l)，主要巖性為淺灰、深灰色白云巖、泥云巖，頂部夾頁巖、含磷介殼灰?guī)r、磷塊巖。

2)奧陶系。紅花園組(Oh)：灰、深灰色生物碎屑灰?guī)r、鮞粒灰?guī)r；湄潭組(Om)：上部灰黃色長石石英砂巖、砂巖，中部灰色介殼灰?guī)r、石英砂巖，下部灰色頁巖夾生物屑灰?guī)r透鏡體；艾家山組(Oa)：上部黃灰色頁狀泥灰?guī)r、泥質(zhì)灰?guī)r；中部鮞?；?guī)r，下部泥質(zhì)生物屑灰?guī)r；寶塔組(Ob)：灰、深灰色泥質(zhì)網(wǎng)紋狀灰?guī)r，頂部為瘤狀灰?guī)r、泥灰?guī)r；龍馬溪組(O3-S1l)：深灰、黑色筆石頁巖夾介殼泥灰?guī)r。

3)志留系。橋溝組(Sq)：深灰色薄層狀泥灰?guī)r與泥、頁巖不等厚互層；石牛欄組(Ss)：上部生物碎屑灰?guī)r，下部瘤狀泥灰?guī)r夾珊瑚層礁；韓家店組(Sh)：灰、黃綠色頁巖、粉砂巖夾生物碎屑灰?guī)r及長石石英砂巖透鏡體。

4)二疊系。梁山組(Pl)：摩尼一帶上部深灰色、灰色泥質(zhì)粉砂巖，下部炭質(zhì)頁巖、灰白色粘土巖，東部僅發(fā)育泥質(zhì)粉砂巖；棲霞組(Pq)：灰、深灰色燧石結(jié)構(gòu)生物屑灰?guī)r，下部眼球狀灰?guī)r；茅口組(Pm)：淺灰、深灰色燧石條帶生物碎屑灰?guī)r，中上部夾灰黑色硅質(zhì)巖，下部眼球。

5)第四系(Qesl)。殘坡積層：棕紅、黃褐色巖塊、亞砂土、亞粘土。

研究區(qū)土地利用類型大部分為旱地、其次水田和果園(圖3)。研究區(qū)土壤環(huán)境和養(yǎng)分指標背景值普遍高于全國土壤(A層)和成都經(jīng)濟區(qū)背景值(表1)。

圖3 研究區(qū)土地利用圖

表1 研究區(qū)土壤相關(guān)元素地球化學特征參數(shù)統(tǒng)計表[3]

2 空間插值方法與強度分析原理

2.1 反距離加權(quán)法(IDW)

在統(tǒng)計中，通常用一組數(shù)據(jù)的均值代表研究對象的真實值，但由于數(shù)據(jù)組中各數(shù)據(jù)對真實值的貢獻大小不同，均值不能較好地反應真實值，因而增加了權(quán)重的概念，即加權(quán)平均法。該種方法增加了各數(shù)據(jù)對真實值的貢獻大小(權(quán)重大小)后求平均值。

對于空間數(shù)據(jù)，離插值點越近的空間樣本點數(shù)據(jù)越接近插值點，表現(xiàn)為離插值點越近的樣本點所占權(quán)重越高。權(quán)重受到距離的影響，且距離越近權(quán)重越高，因而采用加權(quán)平均法對未知點插值，被稱為反距離加權(quán)平均法。實現(xiàn)方法為以插值點與樣本點間的距離為權(quán)重計算插值點S0處的加權(quán)平均值，見式(1)。

(1)

ω(si)=‖si-s0‖-p

(2)

(3)

在對研究區(qū)域內(nèi)進行插值時，各樣本點權(quán)重ω(si)選擇式(2)進行計算，‖si-s0‖表示歐幾里得距離，p為反距離加權(quán)冪，本文使用p為2。與克里金插值相比，反距離插值只考慮了樣本和插值點的距離，而忽略了樣本點的空間結(jié)構(gòu)。

2.2 克里金插值法(Kriging)

空間數(shù)據(jù)往往存在自相關(guān)關(guān)系，且越相鄰的越相似，含量會隨著空間距離的變化呈現(xiàn)規(guī)律的變化。如成礦元素和伴生元素在原生暈或次生暈中存在規(guī)律的變化，靠近礦體元素含量高，遠離礦體元素含量低；同時受到大地構(gòu)造、地質(zhì)背景等影響，區(qū)域上可能呈帶狀分布，而非距離越近含量越高，因而只考慮距離的反距離加權(quán)法不能較好地反應真實情況?？死锝鸩逯捣ú粌H考慮了空間數(shù)據(jù)受距離的影響，也考慮了空間數(shù)據(jù)間的相關(guān)性結(jié)構(gòu)，更真實地反應空間數(shù)據(jù)的分布情況?？死锝鸩逯捣ㄍㄟ^構(gòu)建半變異函數(shù)模型來實現(xiàn)，方法如下：

(4)

圖4 半變異函數(shù)擬合模型

2.3 強度分析原理

在GIS領(lǐng)域，用土地利用強度分析法(LULCC)分析各年度土地利用覆蓋遙感數(shù)據(jù)的年度、土地類型和過度水平變化強度情況，是基于通過各年度土地利用柵格數(shù)據(jù)間交叉表，計算各年度間交叉表數(shù)據(jù)的年度水平、類型水平和過渡水平的土地轉(zhuǎn)變強度來實現(xiàn)[4-6]。筆者將強度分析做了一些改動(修改了R語言包“intensity.analysis”中的部分函數(shù)算法)，基于各插值方法評價后的評價單元柵格數(shù)據(jù)間的交叉表(交叉表反應了兩種插值方法由第一種插值方法轉(zhuǎn)變?yōu)榈诙N插值方法后，前一種插值評價結(jié)果等級轉(zhuǎn)變?yōu)楹笠环N插值評價結(jié)果等級單元個數(shù))，計算交叉表數(shù)據(jù)的插值方法間、評價等級間和等級增量來源情況的轉(zhuǎn)變強度。式(5)為計算t類插值方法間的評價結(jié)果的差異強度，用St表示，其中t表示插值方法組tn-m(前者n類方法轉(zhuǎn)變到后者m類方法)，Ctij代表t類插值方法間(例如，t1分層克里金插值法轉(zhuǎn)變到t2反距離加權(quán)平均插值法)交叉表中，由評價等級i轉(zhuǎn)變?yōu)榈燃塲的評價單元(柵格)個數(shù)。式(6)為計算所有插值方法間的評價結(jié)果差異的平均強度，用U表示；式(7)為計算t類插值方法間中后者插值方法(tm)j類評價等級增加的平均比例，用Gtj表示；式(8)為計算t類插值方法間前者插值方法(tn)i類評價等級減少的平均比例，用Lti表示；式(9)為計算t類插值方法間i類評價等級增加情況，是來自于非i類評價等級轉(zhuǎn)變的結(jié)果，用Rtin表示 ；式(10)為計算t類插值方法間n類評價等級，由其他等級轉(zhuǎn)變而來的平均增加情況，用Wtn表示。

(5)

(6)

(7)

(8)

(9)

(10)

3 空間插值與土地質(zhì)量地球化學評價

為了更好地反應評價單元的真實情況，應盡量保證評價單元內(nèi)有樣品數(shù)據(jù)，本研究調(diào)查的最小采樣單元為250 m×250 m格子，因而將其作為本次的評價單元。評價單元共1 908個，其中“有樣”單元1 794個，占比94.03%；“無樣”單元114個，占比5.97%(圖5)。

圖5 研究區(qū)評價單元分布

圖6 半變異函數(shù)值空間分布圖

圖7 各插值方法土壤N元素分級評價圖

圖8 各插值方法土壤P元素分級評價圖

圖9 各插值方法土壤K元素分級評價圖

土壤質(zhì)量地球化學綜合評價涉及指標包括：養(yǎng)分元素指標N、P、K2O，環(huán)境元素指標As、 Cd、 Cr、 Cu、 Hg、 Ni、 Pb、 Zn以及土壤pH。現(xiàn)主流評價軟件均采用原值賦原評價單元的方式，對“無樣單元”采取插值或賦值的方式。為了對比不同插值方法，增加了單一插值的方式，此外使用主流方法時，對2個及以上樣品數(shù)據(jù)的“有樣”單元，采用了“局部反距離加權(quán)”的方式進行賦值。共涉及8種賦值方式，見method1-8。

method1、method2分別為按地層屬性和按土地利用類型分層克里金插值法，原理為：除了考慮了各元素的半變異函數(shù)模型差異與各向異性外，還兼顧了不同分類區(qū)域存在不同的半變異函數(shù)(協(xié)方差結(jié)構(gòu))；method3為多變量協(xié)同克里金法，原理為：元素不僅自相關(guān)，也與其他變量存在協(xié)相關(guān)，因而該方法構(gòu)建了自相關(guān)和協(xié)相關(guān)半變異函數(shù)模型，但需要設定這些變量具有相同的半變異函數(shù)模型和正定的偏基臺值矩陣，所以采用各元素最優(yōu)半變異函數(shù)模型參數(shù)的中位數(shù)，對評價單元網(wǎng)格進行插值[10]。

Method4為反距離加權(quán)平均法(權(quán)重為式(2))。method5至method8為主流賦值方法，其對“有樣”單元均采用局部加權(quán)平均(權(quán)重為式(3))賦值，而對“無樣”單元插值分別采用“地層屬性分層克里金插值”、“土地利用類型分層克里金插值”、“多變量協(xié)同克里金插值”和“反距離加權(quán)平均法(權(quán)重式(2))插值”。

克里金法插值前，需構(gòu)建半變異函數(shù)模型，其半變異函數(shù)值空間分布圖(圖6)可知： As半變異函數(shù)值在NW-SE方向上含量低于0.4，呈NW-SE向條帶分布，各方向上存在明顯的差異性；其余元素指標半變異函數(shù)值在各方向變化基本一致。因此，為了獲得更準確地評價元素指標的克里金半變異函數(shù)模型，As較其余元素指標需要多考慮方向差異性。最終獲得As在135°方向變程為3 310，在45°方向的變程為3 310×0.7=2 317，其余元素指標半變異函數(shù)模型參數(shù)見表2。method1、method2、method5、method6采用了不同的半變異函數(shù)模型參數(shù)并考慮了As元素方向差異性(模型參數(shù)見表2)；method3、method7使用的半變異函數(shù)模型參數(shù)為各元素最優(yōu)偏基臺和變程的中位數(shù)，分別為0.05和1 584，計算模型為球狀模型“Sph”。

表2 半變異函數(shù)模型最優(yōu)參數(shù)表

根據(jù)《土壤環(huán)境質(zhì)量 農(nóng)用地土壤污染風險管控標準(試行)》(GB15618-2018)和《土地質(zhì)量地球化學評價規(guī)范》(DZ/T 0295-2016)規(guī)范，將最優(yōu)模型插值后的數(shù)據(jù)，進行單指標養(yǎng)分元素N、P、K和土壤質(zhì)量地球化學綜合評價。評價結(jié)果見圖7～圖10，養(yǎng)分單指標和綜合評價結(jié)果的空間分布主要受地質(zhì)背景控制，各插值方法的評價結(jié)果，具有高度一致性。K元素分級結(jié)果很單一，顯示研究區(qū)大面積富K元素。N和P元素分級結(jié)果差異最為明顯，綜合評價結(jié)果次之，表現(xiàn)出N元素豐富、P元素較缺乏和綜合評價的優(yōu)質(zhì)土壤區(qū)域與志留系高橋組分布相一致；局部少數(shù)綜合評價劣等土壤受人類活動影響，分布在研究區(qū)西南處摩尼鎮(zhèn)附近。同時，表現(xiàn)出各指標不同插值法結(jié)果，局部存在明顯差異，集中體現(xiàn)在研究區(qū)的中北部。雖然“無樣”單元比率不足6%，但不同插值方法造成的評價結(jié)果面積的誤差可能很大，最終影響評價結(jié)果的可信度以及政策制定。為了更好的進一步探討不同插值方法間的差異，使用差異度分析方法。

圖10 各插值方法土壤質(zhì)量地球化學綜合評價圖

4 差異度分析

采用R語言“intensity.analysis”代碼包，對各插值法對照組依次采用差異度分析、等級轉(zhuǎn)變強度分析和增量強度來源分析。method1至method8在圖中分別用Ⅰ至Ⅷ的羅馬數(shù)值表示。

4.1 差異度分析

根據(jù)式(5)、式(6)得圖11～圖14。

圖11 氮元素各插值方法差異度圖

圖12 磷元素各插值方法差異度圖

圖13 鉀元素各插值方法差異度圖

圖14 土壤綜合各插值方法差異度圖

總體上鉀元素不同插值法間差異率最低，普遍低于2%和圖8吻合；土壤質(zhì)量地球化學綜合評價結(jié)果各插值法間的差異率較低，普遍低于5%；而單指標氮和磷元素除了method5～method8間，其余方法間差異率普遍高于10%，同圖6和圖7直觀感受一致。

4種指標method1、method2、method3兩兩間差異值都在平均值附近，而method3與method4間差異度遠高于平均值。這說明：克里金法和反距離加權(quán)平均法插值的評價結(jié)果差異很大，而克里金插值法的不同類型方法差異相對較小，但也存在明顯差異。

method1與method5、method2與method6、method3與method7相比較的差異度遠高于平均差異強度，存在強烈的差異。method4與method8相比較的差異度與平均差異度相當，但也存在明顯的差異。這說明：單一用地層分層克里金法、土地利用分層克里金法或多變量協(xié)同克里金法，分別與不同克里金插值法+反距離加權(quán)平均法的差異度遠高于平均差異強度。

4種指標的method5與method6、method6與method7、method7與method8的相比較的差異率不足1%，差異度遠低于平均差異度，這說明“有樣”單元賦反距離加權(quán)平均值與“無樣”單元賦不同方法的插值結(jié)果，都能有效降低方法插值間造成的誤差。

4.2 土壤質(zhì)量地球化學綜合等級轉(zhuǎn)變強度分析

4種指標評價結(jié)果中，最受關(guān)注的為綜合評價指標，為了探討其等級轉(zhuǎn)變方向，因而進一步討論土壤質(zhì)量地球化學綜合等級轉(zhuǎn)變強度。根據(jù)式(7)和式(8)所得各方法組土壤質(zhì)量地球化學綜合等級改變強度圖(圖15)。圖15中紅色條形表示前一種插值方法土地質(zhì)量地球化學綜合等級減少強度，綠色條形表示后一種插值方法土地質(zhì)量地球化學綜合等級增加強度；藍色虛線表示前一種插值方法轉(zhuǎn)變?yōu)楹笠环N插值方法后所有等級增加的平均強度，虛線左側(cè)dormant，表示轉(zhuǎn)變強度不明顯的土地質(zhì)量地球化學綜合等級，虛線右側(cè)active表示轉(zhuǎn)變強度明顯的土地質(zhì)量地球化學綜合等級。

圖15 各方法組土壤質(zhì)量地球化學綜合等級改變強度圖

圖16 各方法組各土壤質(zhì)量地球化學綜合等級增加來源強度圖

圖17 各方法組各土壤質(zhì)量地球化學綜合等級增加來源強度圖

由圖14可知,method5與method6、method6與method7、method7與method8相比較的土壤質(zhì)量地球化學綜合等級轉(zhuǎn)變強度很低，百分比不到15%；method1與method2、method2與method3相比較的等級轉(zhuǎn)變強度最高達到40%；method3～method4、method1～method5、method2～ method6、method3～ method7、method4～ method8的等級轉(zhuǎn)變強度高達80%。

同時，反應出土壤質(zhì)量地球化學綜合評價的優(yōu)質(zhì)和良好等級在兩方法對比的轉(zhuǎn)變強度很大，都表現(xiàn)出了后者方法在優(yōu)質(zhì)和良好的等級較前者方法得到較強的增加。說明：插值方法的選擇會較強地影響到土壤質(zhì)量地球化學綜合評價中的優(yōu)質(zhì)和良好土壤評價。

4.3 增量強度來源分析

根據(jù)式(9)和式(10)所得各方法組各土壤質(zhì)量地球化學綜合等級增量來源強度圖(圖16)。圖16中每列表示增加的土壤質(zhì)量地球化學綜合等級，從左往右依次為優(yōu)質(zhì)、良好、中等和劣等，每行為相同方法組，紅線表示當前等級增加的平均強度，紅線上方target表示當前等級明顯增加的來源等級，紅線下方avoid表示當前等級不明顯增加的來源等級。

method1與method5、method2與method6、method3與method7、method4與method8相比較，后者方法的優(yōu)質(zhì)增加明顯來自于前者方法的良好，優(yōu)質(zhì)土壤增加強度較高達到10%～50%間，呈現(xiàn)了后者方法優(yōu)質(zhì)等級評價單元數(shù)量明顯增加。說明前者插值方法部分評價單元數(shù)據(jù)被平均化，經(jīng)重新加權(quán)平均賦值后，后者方法表現(xiàn)為優(yōu)質(zhì)土地增加。

Method1與method5相比較，后者方法劣等主要來源于前者的優(yōu)級，也同樣說明前者插值方法部分評價單元數(shù)據(jù)被平均化，經(jīng)重新加權(quán)平均賦值后，后者方法表現(xiàn)為劣等土壤增加。

Method5與method6、method6與method7、method7與method8相比較的等級間轉(zhuǎn)變強度很弱，增加強度普遍在2.5%以內(nèi)，且基本與增加強度平均值相當，說明先插值再對有樣品數(shù)據(jù)單元重新加權(quán)平均賦值后，各方法間土壤質(zhì)量地球化學綜合等級接近一致。

5 結(jié)論

1)幾種插值后的評價結(jié)果都能反映其空間分布受到地質(zhì)背景控制的影響。

2)通過借鑒土地利用強度分析法(LULCC)得知：①差異度分析,“有樣”單元賦反距離加權(quán)平均值與“無樣”單元賦不同方法的插值結(jié)果，都能有效降低方法插值間造成的誤差；②等級轉(zhuǎn)變強度分析,方法間差異度較大的原因在于優(yōu)質(zhì)、良好和劣等間差異度較大，其等級內(nèi)差異度普遍在20%～80%間，插值法+反距離加權(quán)平均結(jié)合法其優(yōu)良和劣等的增加強度明顯，單一插值法其優(yōu)良減少強度明顯；③等級增量來源強度分析,插值法+反距離加權(quán)平均結(jié)合法其優(yōu)質(zhì)、劣等的增加量，主要分別源于單一插值法的良好和優(yōu)質(zhì)；總體呈現(xiàn)單一插值法靠中間的等級減少向插值法+反距離加權(quán)平均法兩端等級增加的方向轉(zhuǎn)變的趨勢。

3)“有樣”單元賦反距離加權(quán)平均值與“無樣”單元賦不同方法的插值結(jié)果相結(jié)合，能避免因空間插值計算范圍過大，促使高值或低值數(shù)據(jù)被平均后，導致較小的評價單元評價等級降低或升高的情況發(fā)生。

4)筆者只評價了250 m×250 m的柵格單元，有必要進一步在不同尺度柵格(Raster)和多邊形(Polgon)圖斑的評價單元上探討，以便使不同尺度調(diào)查評價結(jié)果更趨于真實地反映客觀情況。

致謝

感謝單位同事提供數(shù)據(jù)支持。感謝北京師范大學，韓春萌博士，對本人ARCGIS的使用提供幫助。