基于半變異函數(shù)的土壤有機(jī)碳空間預(yù)測(cè)方法

苗正紅，邱中軍，畢 強(qiáng)，楊清臣

（1.吉林省水利水電勘測(cè)設(shè)計(jì)研究院，吉林 長(zhǎng)春 130021）

土壤圈是一個(gè)不均衡系統(tǒng)，元素的遷移、分散和累積是受地域分異規(guī)律所制約的。土壤有機(jī)碳近年來(lái)成為人們關(guān)注的熱點(diǎn)。目前，國(guó)內(nèi)外有關(guān)土壤有機(jī)碳空間預(yù)測(cè)的方法主要有土壤類(lèi)型法、植被類(lèi)型法/生態(tài)系統(tǒng)類(lèi)型法、反距離加權(quán)方法等[1-4]。由于土壤分布的空間變異性和各區(qū)域相關(guān)因素的差異性，使得這些方法在研究中受到不同程度的限制，其準(zhǔn)確性和空間模擬效果不盡理想，且應(yīng)用性較差。本文基于半變異函數(shù)，對(duì)三江平原土壤有機(jī)碳數(shù)據(jù)進(jìn)行空間預(yù)測(cè)，并將其結(jié)果與反距離加權(quán)法空間預(yù)測(cè)的結(jié)果進(jìn)行比較，對(duì)比RMSE和空間模擬效果，從而提出一種快速、準(zhǔn)確、應(yīng)用性較好的土壤有機(jī)碳空間預(yù)測(cè)方法。

1 材料與方法

1.1 研究區(qū)概況

三江平原包括完達(dá)山以北的松花江、黑龍江和烏蘇里江沖積形成的低平原和完達(dá)山以南的烏蘇里江及其支流與興凱湖形成的沖積湖積平原，總面積108 900 km2，平均海拔50～60 m，地勢(shì)由西南向東北緩緩傾斜。

1.2 數(shù)據(jù)預(yù)處理

在三江平原選取水田、旱田、林地、草地、濕地5種類(lèi)型土地覆被類(lèi)型，土地利用/覆蓋數(shù)據(jù)由2010～2011年的Landsat TM遙感數(shù)據(jù)經(jīng)預(yù)處理、幾何糾正、圖像鑲嵌拼接和裁剪，并經(jīng)人工目視解譯得到，已進(jìn)行廣泛的地面調(diào)查、咨詢(xún)、驗(yàn)證和精度評(píng)價(jià)。2010年在三江平原調(diào)查了367處樣地（圖1），在每個(gè)采樣點(diǎn)周?chē)?個(gè)點(diǎn)，混合土樣，四分法取樣，其中濕地采樣方法是在每個(gè)采樣點(diǎn)選擇3 塊4 m ×4 m的樣地，采樣前先去除地表未腐解的凋落物，采用多點(diǎn)混合采樣法采集0～30 cm深度的土壤樣品，裝入無(wú)菌袋中。采集的土樣迅速帶回實(shí)驗(yàn)室，剔除可見(jiàn)的動(dòng)、植物殘?bào)w和石塊，用四分法取出適量土壤樣品用于測(cè)定土壤有機(jī)碳。土壤容重采用環(huán)刀法（100 cm3）采集0～30 cm的土壤樣品，所采集的土壤樣品自然風(fēng)干后，以四分法取樣。用于測(cè)定有機(jī)碳的土樣，風(fēng)干后挑去根系和大于2 mm的礫石，利用FW-100粉碎機(jī)粉碎，過(guò)100目篩。本研究采用重鉻酸鉀外加熱氧化法測(cè)定土壤有機(jī)碳含量[5]。

圖1 三江平原2010年土壤采樣點(diǎn)與土地利用空間分布圖（審圖號(hào)：黑S（2010）051號(hào)）

1.3 半變異函數(shù)和反距離加權(quán)方法

1.3.1 半變異函數(shù)

半變異函數(shù)是一個(gè)描述土壤特性空間變異結(jié)構(gòu)的函數(shù)，通過(guò)系數(shù)R2來(lái)決定并綜合考慮RSS（殘差）、塊金值和有效距離確定最優(yōu)的半變異函數(shù)理論模型：

式中，r（h）是樣本距為h的半方差；h是樣本距（變程lag）；N（h）是間距為h的樣本對(duì)的總個(gè)數(shù)；Z（xi）是樣本z在位置xi的實(shí)測(cè)值，Z（xi+h）是與xi距離為h處的樣本值。

Kriging方法根據(jù)半方差函數(shù)分析提供的空間自相關(guān)程度的信息進(jìn)行插值，將未測(cè)定處x0的估計(jì)值Z'（x0）假設(shè)為已知觀測(cè)值的線性和[6]：

式中，λi為與測(cè)點(diǎn)有關(guān)的加權(quán)系數(shù)。基于半方差函數(shù)，Kriging插值可根據(jù)無(wú)偏估計(jì)和方差最小兩項(xiàng)要求來(lái)確定。

1.3.2 反距離加權(quán)方法

反距離加權(quán)法是一種常用的簡(jiǎn)單的空間插值方法，它是基于“地理第一定律”的基本假設(shè)，即兩個(gè)物體相似性隨二者間距離的增大而減少。

設(shè)有n個(gè)點(diǎn)，平面坐標(biāo)為（xi，yi），垂直高度為zi，i=1，2，…，n時(shí)，倒數(shù)距離加權(quán)插值的插值函數(shù)為：

式中，當(dāng)（x，y）≠（xi，yi），i=1，2，…，n時(shí)，f（x，y）（x，y）=（xi，yi），i=1，2，…，n時(shí)，f（x，y）=zj。

反距離加權(quán)插值是以插值點(diǎn)與樣本點(diǎn)間的距離為權(quán)重進(jìn)行加權(quán)平均的，離插值點(diǎn)越近的樣本賦予的權(quán)重越大，該方法簡(jiǎn)單、直觀，且效率高。

1.4 基于半變異函數(shù)空間預(yù)測(cè)的技術(shù)路線

圖2為半變異函數(shù)空間預(yù)測(cè)方法的技術(shù)路線，首先確定土壤有機(jī)碳數(shù)據(jù)是否符合正態(tài)分布，這需要在SPSS中進(jìn)行描述性分析，通過(guò)P_P圖（圖3）和表1可以判斷該數(shù)據(jù)是否符合正態(tài)分析，如果符合，則在GS+中確定數(shù)據(jù)適應(yīng)的模型，然后計(jì)算模型最佳指數(shù)（IGF）是否在最佳范圍內(nèi)。如果為最佳，則按照GS+設(shè)定的參數(shù)在ArcGIS中進(jìn)行空間插值預(yù)測(cè)；否則重新設(shè)定參數(shù)。

本文以三江平原土壤數(shù)據(jù)為例，利用GS+和Variowin軟件來(lái)計(jì)算半變異函數(shù)及確定最優(yōu)模型[7,8]，模型的IGF在0.000 5～0.005 9之間，表示該模型模擬效果最佳。采用ArcGIS 9.3和Origin 8.0軟件進(jìn)行空間插值和制圖。

圖2 技術(shù)路線

表1 三江平原表層（0～30 cm）土壤有機(jī)碳密度描述性統(tǒng)計(jì)

圖3 三江平原表層土壤有機(jī)碳密度P-P圖

2 結(jié)果分析

2.1 土壤有機(jī)碳密度描述性特征

根據(jù)表1中對(duì)偏度系數(shù)和峰度系數(shù)的計(jì)算，并從圖3對(duì)比中可以看出，2010年表層土壤有機(jī)碳密度（SOCD）不符合正態(tài)分布，但是對(duì)數(shù)轉(zhuǎn)換后符合正態(tài)分布，所以可以進(jìn)行空間插值。

2.2 半變異函數(shù)參數(shù)確定

2010年表層土壤有機(jī)碳密度的半變異函數(shù)是高斯模型，從表2中可以看出，IGF值為0.004 47，屬于最佳擬合（圖4）。2010年的塊金效應(yīng)為15%（<25%），表示其具有強(qiáng)烈的空間相關(guān)性，主要是由于土壤母質(zhì)、地形、土壤類(lèi)型等非人為的結(jié)構(gòu)性因素引起的變異。

圖4 半變異函數(shù)圖

表2 表層（0～30 cm）土壤有機(jī)碳密度半方差模型參數(shù)

2.3 兩種方法空間預(yù)測(cè)結(jié)果對(duì)比與驗(yàn)證

2.3.1 空間預(yù)測(cè)分布圖對(duì)比

從空間上看，半變異函數(shù)的預(yù)測(cè)結(jié)果空間分布過(guò)度平滑，無(wú)異常值，整體趨勢(shì)特征明顯；而反距離加權(quán)方法預(yù)測(cè)結(jié)果空間分布過(guò)度生硬，整體趨勢(shì)不夠明顯。說(shuō)明基于半變異函數(shù)的空間預(yù)測(cè)結(jié)果要優(yōu)于反距離加權(quán)法。

圖5 兩種方法空間預(yù)測(cè)結(jié)果對(duì)比

2.3.2 預(yù)測(cè)精度對(duì)比

本文利用研究區(qū)未參與空間預(yù)測(cè)的38個(gè)野外實(shí)測(cè)數(shù)據(jù)作為驗(yàn)證點(diǎn)數(shù)據(jù)，通過(guò)驗(yàn)證點(diǎn)的實(shí)測(cè)值和預(yù)測(cè)值計(jì)算不同方法的RMSE和平均誤差，結(jié)果如表3。從表3可以看出，半變異函數(shù)的空間預(yù)測(cè)方法的RMSE和平均誤差均比反距離加權(quán)法的結(jié)果更接近0，這表明半變異函數(shù)的土壤有機(jī)碳空間預(yù)測(cè)方法預(yù)測(cè)精度更高。

表3 不同預(yù)測(cè)方法預(yù)測(cè)精度比較

為了直觀反映出預(yù)測(cè)精度的對(duì)比情況，本文以反距離加權(quán)法的RMSE為基準(zhǔn)，計(jì)算得到基于半變異函數(shù)的預(yù)測(cè)精度的相對(duì)提高值，結(jié)果提高14%。總體而言，半變異函數(shù)的空間預(yù)測(cè)精度要高于反距離加權(quán)預(yù)測(cè)精度，半變異函數(shù)更適合土壤有機(jī)碳空間預(yù)測(cè)。

3 結(jié) 語(yǔ)

本文基于半變異函數(shù)原理，以三江平原土壤數(shù)據(jù)為例，對(duì)土壤有機(jī)碳進(jìn)行空間預(yù)測(cè)，其預(yù)測(cè)結(jié)果與反距離加權(quán)方法的結(jié)果進(jìn)行對(duì)比。結(jié)果表明，基于半變異函數(shù)的預(yù)測(cè)結(jié)果優(yōu)于反距離加權(quán)法，且基于半變異函數(shù)的空間預(yù)測(cè)法不受地域限制，其應(yīng)用性較強(qiáng)，方法便捷，為今后其他土壤數(shù)據(jù)快速準(zhǔn)確地空間預(yù)測(cè)提供基礎(chǔ)。

[1]王紹強(qiáng),劉紀(jì)遠(yuǎn),于貴瑞.中國(guó)陸地土壤有機(jī)碳蓄積量估算誤差分析[J].應(yīng)用生態(tài)學(xué)報(bào),2003,14(5):797-802

[2]Post W M,Emanuel W R,Zike P J,et al.Soil Carbon Pools and World Life Zones[J].Nature,1982,298: 56-159

[3]楊元合.青藏高原高寒草地生態(tài)系統(tǒng)碳氮儲(chǔ)量[D].北京:北京大學(xué),2008

[4]Levesque J,King D J.Airborne Digital Camera Image Semivariance for Evaluation of Forest Structural Damage at an Acid Mine Site[J].Remote Sensing of Environment,1999,68:112-124

[5]Tiessen H,Moir J O.Total and Organic Carbon[J].Soil Sampling and Methods of Analysis,1993:187-199

[6]王政權(quán).地統(tǒng)計(jì)學(xué)及在生態(tài)學(xué)中的應(yīng)用[M].北京:科學(xué)出版社,1999

[7]舒彥軍,張立亭.半變異函數(shù)球狀模型2次套合結(jié)構(gòu)優(yōu)化算法[J].地理空間信息,2012,10(2):61-63

[8]Pannatier Y.Variowin:Software for Spatial Data Analysis in 2D[M].New York:Springer-Verlag,1996