多種空間插值法在土地質量地球化學評價中對比分析
——以四川省敘永縣麻城鎮和摩尼鎮為例

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

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

0 引言

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

1 研究區概況

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

圖1 研究區地理位置圖[3]

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

圖2 研究區地質圖[3]

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

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

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

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

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

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

圖3 研究區土地利用圖

表1 研究區土壤相關元素地球化學特征參數統計表[3]

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

2.1 反距離加權法(IDW)

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

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

(1)

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

(2)

(3)

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

2.2 克里金插值法(Kriging)

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

(4)

圖4 半變異函數擬合模型

2.3 強度分析原理

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

(5)

(6)

(7)

(8)

(9)

(10)

3 空間插值與土地質量地球化學評價

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

圖5 研究區評價單元分布

圖6 半變異函數值空間分布圖

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

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

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

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

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

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

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

表2 半變異函數模型最優參數表

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

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

4 差異度分析

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

4.1 差異度分析

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

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

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

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

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

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

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

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

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

4.2 土壤質量地球化學綜合等級轉變強度分析

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

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

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

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

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

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

4.3 增量強度來源分析

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

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

Method1與method5相比較，后者方法劣等主要來源于前者的優級，也同樣說明前者插值方法部分評價單元數據被平均化，經重新加權平均賦值后，后者方法表現為劣等土壤增加。

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

5 結論

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

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

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

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

致謝

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