基于回歸克里格的河北省土壤有機質(zhì)時空變化特征研究

李 可，楊 勇，劉亞軍，譚秀麗，楊 雪，李恩光，吳若景

（華中農(nóng)業(yè)大學(xué)資源與環(huán)境學(xué)院，農(nóng)業(yè)農(nóng)村部長江中下游耕地保育重點實驗室，湖北 武漢 430070）

土壤是人類進行各項生產(chǎn)和賴以生存的基地，土壤質(zhì)量變化與質(zhì)量高低關(guān)系到人類的生存狀況，而土壤有機質(zhì)（SOM）是構(gòu)成土壤有機無機復(fù)合膠體的核心物質(zhì)，是土壤的重要組成部分，對土壤結(jié)構(gòu)、物理性狀、化學(xué)性狀和生態(tài)過程有著重要的影響與作用［1-2］。所以，研究SOM的時空變化特征，對協(xié)調(diào)耕地土壤肥力，提高土壤質(zhì)量水平，揭示土壤有機質(zhì)空間分布格局，了解土壤質(zhì)量動態(tài)變化，優(yōu)化農(nóng)業(yè)種植結(jié)構(gòu)和施肥管理具有科學(xué)指導(dǎo)作用［3-4］。

近年來，由于GIS技術(shù)和地統(tǒng)計學(xué)的快速發(fā)展，學(xué)者們開始利用兩者相結(jié)合來研究各個尺度上SOM時空變化特性，為測土培肥等技術(shù)實施提供理論支持。如田塊尺度上，李國良等［5］探究了廣東荔枝園土壤養(yǎng)分肥力時空變化特征，提出該地區(qū)土壤養(yǎng)分肥力發(fā)生了很大變化，多種肥力屬性有明顯下降趨勢，但有機質(zhì)變化不大；縣市尺度上，王少賀等［6］對臺安縣SOM的時空變化特征的研究表明：25年來，該縣耕地土壤有機質(zhì)平均含量下降了2.65 g/kg，整體已呈現(xiàn)出下降的趨勢，這與人們長期的“重種輕養(yǎng)、只用不養(yǎng)”有關(guān)；鮑思屹等［7］研究了建甌市近40年來竹林土壤有機質(zhì)的變化特征，得出建甌市毛竹林土壤有機質(zhì)Ⅰ級與Ⅲ級土壤面積分別增加了4.68%和17.3%，而Ⅱ與Ⅳ級土壤面積分別減少了17.4%和4.63%，人類活動對竹林土壤有機質(zhì)的影響程度有所增強；省級尺度上，趙明松等［8］以江蘇省為研究對象，對比分析了近26年全省土壤有機質(zhì)的時空變異特征，結(jié)果表明：SOM含量空間分布呈現(xiàn)出北增南減，沿江平原增，寧鎮(zhèn)丘陵減，濱海平原基本持平的空間格局，增加幅度由北向南逐漸減??；全國尺度上，任意等［9］對我國東北、西北、華北、西南、中南和華東不同地區(qū)土壤養(yǎng)分的差異及變化趨勢進行了研究，結(jié)果指出在1986～1997年中南、華北區(qū)有顯著的上升趨勢，華東、東北區(qū)有下降趨勢，西北、西南區(qū)變化平穩(wěn)，1998～2006年各個區(qū)域的土壤有機質(zhì)變化基本平穩(wěn)。以上研究均表明，SOM在各空間尺度上存在明顯的時空分異，這對分析土壤肥力對人類活動的響應(yīng)、指導(dǎo)農(nóng)業(yè)生產(chǎn)、實現(xiàn)農(nóng)業(yè)生態(tài)的可持續(xù)發(fā)展有重要意義。

從地理條件上講，河北省地形地貌復(fù)雜，是全國唯一兼有高原、山地、丘陵、平原、湖泊和海濱的省份，土壤類型較多。從經(jīng)濟人口條件上講，河北省總?cè)丝跀?shù)全國排第6，約為7 185.4萬（第六次全國人口普查數(shù)據(jù)）。而且河北省是農(nóng)業(yè)大省，全省34.60%的面積為耕地，其糧食產(chǎn)量在全國位居前列。因此，本文選取了河北省作為研究區(qū)域，并采用GIS空間分析、地統(tǒng)計分析法以及回歸克里格空間插值方法，對1980～2010年間河北省土壤有機質(zhì)的時空變化進行了研究，為預(yù)測糧食產(chǎn)量、優(yōu)化施肥管理、保護生態(tài)平衡和耕地的可持續(xù)發(fā)展等提供了科學(xué)依據(jù)。

1 材料與方法

1.1 研究區(qū)概況

河北省簡稱“冀”，位于東經(jīng)113°27′～119°50′，北緯 36°05′～42°40′，地處華北地區(qū)，北靠燕山，南望黃海，西依太行，東臨渤海，面積約18.8萬km2。河北省地勢自西北向東南傾斜，東部有太行山脈，全省依據(jù)地形地貌具體可分為河北平原區(qū)、壩上高原區(qū)、燕山和太行山地區(qū)3大部分。全省最高海拔2 855 m，平均海拔551 m，大部分地區(qū)還是屬于平原丘陵?？偟膩碚f，河北省是全國唯一兼有高原、山地、丘陵、平原、湖泊和海濱的省份，其土壤類型按發(fā)生分類共分為21種土類、55種亞類，整體呈帶狀分布，變化較大，極具研究意義。

1.2 數(shù)據(jù)來源

1978～1980年（以下均表示為1980年）河北省表層土壤數(shù)據(jù)來自全國第二次土壤普查時期，《河北土種志》［10］中記錄的334個典型土壤剖面，每個剖面的坐標(biāo)位置通過《土種志》對典型剖面的地理位置描述與高清遙感影像結(jié)合得出。2009～2011年（以下均表示為2010年）表層土壤數(shù)據(jù)為國家科技基礎(chǔ)性工作專項“我國土系調(diào)查與《中國土系質(zhì)》［11］編制”中河北省162個典型土壤剖面數(shù)據(jù)。兩個時期具體樣點分布如圖1所示，不同的土地利用類型上均有樣點分布，1980年時耕地、林地、草地、水域、未利用地與城鄉(xiāng)、工礦、居民用地6類土地利用類型上分別分布了195、30、31、10、5、63個土壤樣點，2010年時則分別分布了68、41、35、6、4、8個土壤樣點，其中土系調(diào)查72個典型土壤剖面均在第二次土壤普查時期的典型剖面附近采集。

本文采用的輔助變量為海拔高程、曲率、歸一化植被指數(shù)、坡度、土地利用類型和土壤類型。因此，除了土壤剖面數(shù)據(jù)以外，還收集了研究區(qū)域1 km分辨率的數(shù)字高程圖（DEM）、2010年植被指數(shù)（NDVI）空間分布圖、1980和2010年土地利用數(shù)據(jù)（由Landsat TM/ETM遙感圖像人工目視解譯生成）等環(huán)境因子數(shù)據(jù)，來源于中國科學(xué)院資源環(huán)境科學(xué)數(shù)據(jù)中心（http：//www.resdc.cn）。河北省1∶100萬土壤類型圖為ArcGIS數(shù)字化數(shù)據(jù)，來源于國家科技基礎(chǔ)條件平臺—國家地球系統(tǒng)數(shù)據(jù)共享服務(wù)平臺-土壤科學(xué)數(shù)據(jù)中心（http：//soil.geodata.cn）。坡度和曲率數(shù)據(jù)在ArcGIS 10.2軟件中由DEM數(shù)據(jù)提取而來。結(jié)果分析時涉及統(tǒng)計數(shù)據(jù)來源于1995和2015年《河北農(nóng)村統(tǒng)計年鑒》、1985和2011年《河北經(jīng)濟年鑒》，涉及部分土地利用數(shù)據(jù)來自上述土地利用遙感數(shù)據(jù)計算得出的結(jié)果。

圖1 研究區(qū)域土壤樣點分布圖

1.3 研究方法

1.3.1 回歸克里格插值（Regression Kriging，RK）

回歸克里格的基本步驟原理：首先，對輔助環(huán)境變量和目標(biāo)變量之間的回歸關(guān)系進行計算；然后，對空間樣點上的趨勢項進行分離，用實際值減去趨勢項值得到殘差，對殘差進行線性無偏估計處理；最后，將回歸分析的趨勢項同殘差插值結(jié)果相加，從而得到目標(biāo)變量的空間估測值［12-13］。假設(shè)目標(biāo)變量為 z（s1），z（s2），…，z（sn），si=（xi，yi），i=1，2，…，n，si為實測點的坐標(biāo)值，n為實測點的數(shù)目。則預(yù)測點s0的值z（s0）可表示為：

式中，趨勢項m（s0）通常用一元或多元線性回歸進行擬合，殘差項ε（s0）用普通克里格插值。由上式可以看出，回歸克里格方法能夠綜合多個確定性的輔助環(huán)境變量（如曲率、坡度、高程等）對目標(biāo)變量進行解釋預(yù)測［14］，在樣點較稀疏的區(qū)域得到精度較高的預(yù)測結(jié)果［15］。鑒于本研究的區(qū)域面積較大，樣點數(shù)量相對較少，因此，選用此方法對土壤養(yǎng)分進行空間插值。

1.3.2 啞變量引入

土地利用和土壤類型分布圖均為范疇型變量，與DEM、NDVI等連續(xù)型變量不同的是，其Value值分別為標(biāo)識代碼和實際數(shù)值，兩者無法同時直接用于進行回歸分析，因此引入啞變量值［16］。啞變量為虛擬變量，常用于處理定性因子或分類變量，一般取值為0或1［17］。啞變量的定義為對于等級（定性）數(shù)據(jù)α，用變量 （α，β）表示為：

這種方法叫做定性因子二值化展開，變量（α，β）就稱為啞變量。因定性變量為取0或1的數(shù)值向量，便于用數(shù)值方法進行處理［18］。本文將土地利用數(shù)據(jù)按照一級分類分為耕地、林地、草地、水域、未利用地與城鄉(xiāng)、工礦、居民用地6類，土壤類型數(shù)據(jù)按照發(fā)生分類中的土綱分為9類，然后按上述方法進行啞變量換算，符合相應(yīng)類型即為1，不符合即為0，最后將啞變量和其它連續(xù)變量結(jié)合做出回歸分析。另外，本研究是直接利用全國第二次土壤普查圖件河北省土壤圖提取土綱進行啞變量轉(zhuǎn)換，通過地圖套合與2010年土系調(diào)查的土壤樣點對應(yīng)獲得樣點發(fā)生分類名稱，而非將2010年土壤樣點的系統(tǒng)分類的土壤類型直接對照為發(fā)生分類土壤類型。

1.4 分析工具

回歸克里格的計算與插值采用Matlab R2014a軟件和ArcGIS 10.2的地統(tǒng)計模塊來完成；空間數(shù)據(jù)的處理與制圖依靠ArcMap 10.2來完成；基本的統(tǒng)計分析由SPSS 19和Excel 2010軟件完成。

2 結(jié)果與分析

2.1 1980～2010年全省土壤有機質(zhì)的時間變化

表1為第二次全國土壤普查土壤有機質(zhì)分級標(biāo)準(zhǔn)［19］，表2為1980和2010年土壤典型剖面點描述性統(tǒng)計結(jié)果。由表2可知，1980和2010年土壤有機質(zhì)含量范圍分別為0.80～164.50和2.60～208.50 g/kg，均值分別為22.02和18.38 g/kg，屬于III和IV級標(biāo)準(zhǔn)（對比表1），這說明從總體上看，近30年來SOM含量有一定程度的下降。而且，1980和2010年SOM含量變異系數(shù)分別為1.10和1.07，均大于1，屬于高度變異，這表示在該時期內(nèi)河北省的該類土壤有機質(zhì)含量的空間分布差異較大，受人為影響嚴(yán)重。同時，對比兩時期內(nèi)SOM含量的變異系數(shù)可發(fā)現(xiàn)，近30年來土壤有機質(zhì)含量的變異系數(shù)有所減小，這表示隨著耕作制度、施肥量變化、施肥方式和田間管理措施等人為因素的影響，土壤有機質(zhì)含量的空間變異特征表現(xiàn)為逐漸減弱的趨勢。

表1 第二次全國土壤普查土壤有機質(zhì)分級標(biāo)準(zhǔn) （g/kg）

表2 1980和2010年土壤樣點SOM描述性統(tǒng)計數(shù)據(jù)

表3為河北省1980和2010年各地市SOM含量變化結(jié)果。由表3來看，1980年時石家莊、保定、張家口、承德地區(qū)的SOM含量均高于全省平均值，2010年時只有承德地區(qū)SOM含量高于全省平均值。從兩個時期對比變化來看，保定和張家口地區(qū)SOM含量在1980年時明顯高于其他地區(qū)，但隨著近30年來SOM含量變化，其整體呈現(xiàn)降低趨勢，且差值很大，均從II級標(biāo)準(zhǔn)直跌至IV級標(biāo)準(zhǔn)。石家莊、承德、邢臺3地SOM含量同樣呈降低趨勢。其中，石家莊和承德SOM含量降低量適中，邢臺地區(qū)SOM含量降低量較小。其他地市中，除邯鄲地區(qū)SOM含量無明顯變化外，剩余地區(qū)SOM含量均呈增長狀態(tài)，且增量都較小。

造成上述土壤有機質(zhì)時間變化的原因：張家口和承德均位于北部，此地以棕壤和栗鈣土為主，腐殖質(zhì)層較厚，所以兩地區(qū)在1980年時土壤有機質(zhì)較高。而在1990～2008年間，由于大量退耕還林還草導(dǎo)致張家口和承德兩地常用耕地面積分別減少了220 362、79 709 hm2，但近30年中兩地區(qū)總耕地面積增量百分比分別只有0.65%、-0.11%。這表明至2010年，張家口和承德地區(qū)在總耕地面積變化極小的前提下，大量非耕地轉(zhuǎn)變?yōu)楦?，所以SOM呈下降趨勢。同時，張家口地區(qū)耕地壓力指數(shù)增大［20］，至2010年，其耕地肥料使用量雖增加至113.50 kg/hm2，但相對其它市區(qū)來說其施肥量最低，所以土壤養(yǎng)分無法得到及時補充，也會導(dǎo)致SOM含量降低。承德地區(qū)耕地壓力指數(shù)也在增加，加之近些年各類土地利用強度加大，農(nóng)業(yè)技術(shù)及管理水平相對較低，因此土壤有機質(zhì)在大幅度消耗后無法及時補足，SOM含量降低［20-21］。而保定、石家莊、邢臺地區(qū)，1980年SOM含量較高樣點所處的區(qū)域不屬于耕地范圍，所以不會對其進行施肥等處理，隨著土壤養(yǎng)分流失，其SOM含量值同樣也就會降低。而其它地區(qū)可能由于土壤類型、20世紀(jì)80年代以前本省耕地管理問題等的影響［22］，在1980年時土壤有機質(zhì)相較于上述地市偏低，但經(jīng)過近30年施肥、秸稈還田、耕作制度調(diào)整等措施，促進了土壤有機質(zhì)的提升。

表3 1980～2010年河北省各地市土壤有機質(zhì)含量對比

2.2 回歸克里格插值結(jié)果

利用回歸克里格插值方法對兩時期樣點數(shù)據(jù)進行插值。在插值過程中分別隨機選取20%樣點作為驗證點進行精度評價，不參與插值，回歸分析及精度評價結(jié)果見表4。表4中R2、F、P為回歸方程評價指標(biāo)，RMSE為插值精度評價指標(biāo)。由表4可得，1980和2010年SOM的分析P值均小于0.01，屬于極顯著水平。雖然SOM在兩時期內(nèi)的決定系數(shù)和統(tǒng)計量較低，但其最終兩時期SOM含量的結(jié)果精度值RMSE均小于其對應(yīng)均值的1/2（對比表2），所以整體插值精度也屬于誤差接受范圍。特別的是，兩時期中1980年SOM含量的RMSE值更小，可能是由于此時期用于插值的樣點數(shù)目為2010年的2倍，從而使得插值結(jié)果的精度更高。

2.3 1980～2010年全省土壤有機質(zhì)的空間變化

圖2為1980～2010年SOM含量空間分布圖。由圖2可知，1980和2010年SOM含量整體上空間分布均是西北高東南低，并且整體呈現(xiàn)由西北向東南降低的趨勢，至2010年SOM含量東南與西北差異減小。由圖2a可得，1980年時SOM含量達(dá)到I、II級標(biāo)準(zhǔn)的高值區(qū)域主要分布在承德西北部、張家口東北部和南部、保定西北部，V、VI級的低值區(qū)域主要在衡水、邢臺東部、唐山北部。且SOM值大多集中在III級和IV級標(biāo)準(zhǔn)區(qū)間內(nèi)，分別占全省面積的16.35%和33.85%，主要分布在河北省南部平原地區(qū)。由圖2b分析得，2010年SOM含量高值區(qū)域主要在承德、張家口西北部、保定北部等地，而低值區(qū)域面積較小，分布零散。此時，SOM含量同樣集中在III級和IV級標(biāo)準(zhǔn)區(qū)間內(nèi)，分別占全省面積的47.82%和20.73%，主要分布在南部平原地區(qū)、唐山地區(qū)。

圖3為兩時期SOM含量的空間變化分布圖，白色和淺灰色區(qū)域代表SOM含量減少，深灰色和黑色區(qū)域代表其含量增加。從整體上看，SOM含量整體上呈現(xiàn)西北減東南增的趨勢，增幅總體上由西北向東南依次增加；SOM含量增加的面積占全省面積的65.54%，主要分布在南部平原、唐山、秦皇島地區(qū)，其中增幅在0～10 g/kg間占比最高，約占全省面積的36.29%。從局部看，SOM減少區(qū)域主要位于1980年SOM高值分布地區(qū)，SOM增值較高的區(qū)域則對應(yīng)處于1980年SOM低值分布區(qū)域。同時，SOM含量增幅最大（＞20 g/kg）的區(qū)域在承德東部、唐山北部和秦皇島，降幅最大（＜-20 g/kg）的區(qū)域在承德西北部、石家莊和保定西部。造成上述變化的原因是由于降幅較大的地區(qū)在1980年時SOM含量就較高，但多為草地、林地，基本不屬于耕地。其中，位于承德西北部的圍場、豐寧地處河北壩上地區(qū)，至2005年時，此處放牧、樵采過度，壩上草場面積已由20世紀(jì)50年代初的86.7萬hm2減少到現(xiàn)在的22.6萬hm2，草場覆蓋度由90%降低到44%，森林面積也呈明顯減少趨勢，土地沙漠化嚴(yán)重，導(dǎo)致該區(qū)域土壤肥力流失嚴(yán)重，SOM含量大幅度降低［23］。而石家莊和保定西部地處太行山區(qū)，屬于石質(zhì)山區(qū)，由于過度追求開墾種植、砍伐森林，導(dǎo)致森林覆蓋率低，生態(tài)環(huán)境脆弱，土層養(yǎng)分涵養(yǎng)差，SOM含量顯著降低［24］。而原來的低值區(qū)域大多處于耕地范圍，自20世紀(jì)90年代起，河北省開始大力推廣綠肥種植、秸稈還田等措施，目前河北省東南部平原區(qū)域基本實現(xiàn)100%秸稈還田率，使得該區(qū)域土壤有機質(zhì)狀況得到較好改善［21，25］。

表4 1980和2010年河北省SOM回歸分析及插值精度結(jié)果

圖2 1980～2010年SOM空間分布圖

圖3 1980～2010年SOM空間變化分布圖

3 討論

由土壤有機質(zhì)質(zhì)量的演化過程可知，1980～2010 年土壤有機質(zhì)含量的時空變化主要與自然和社會經(jīng)濟發(fā)展有關(guān)，自然過程的影響因素主要涉及氣候、溫度、降水等，但這些因素的變化從近30年來看，總體波動幅度較小，對SOM的累積與變化影響不大，甚至可忽略不計。因此，本文主要側(cè)重分析社會經(jīng)濟因素對土壤有機質(zhì)變化的影響。由李繼明等［26］的研究可知，綠肥與化肥長期配合施用會使得土壤有機質(zhì)、全氮和全磷均有所積累，其積累的量與肥料施用量及有機肥種類有關(guān)。所以，接下來本研究主要分析土壤有機質(zhì)含量的時空變化同肥料施用量間的關(guān)系。

1980～2010年河北省化肥施用量（折純量）由73.10萬t增加至322.86萬t，耕地面積由664.79萬hm2減少至598.89萬hm2，平均用量由109.96 kg/hm2增加至539.10 kg/hm2，其中氮肥、磷肥、鉀肥和復(fù)合肥的用量有不同程度的增加（圖4）。1985～2010年間，全省氮肥和復(fù)合肥用量增加較大，分別由109.79、24.38 kg/hm2增加至255.64、159.63 kg/hm2；磷肥和鉀肥的用量增加較平穩(wěn)。同時，1980～2010年期間，全省糧食（稻谷、小麥、玉米、大豆、薯類）總產(chǎn)量由1 575萬t增加至2 976萬t，糧食平均產(chǎn)量由2.15 t/hm2增加至4.88 t/hm2，這主要得益于種植高產(chǎn)作物及使用更多的化肥、農(nóng)藥、地膜等［20］。同時，由于化肥使用量的大幅度提升，全省耕地土壤肥力質(zhì)量也有一定提升［22］。而對比SOM空間變化分布圖也可發(fā)現(xiàn)，絕大多數(shù)SOM含量增長的區(qū)域?qū)儆诟胤秶?，與上述分析結(jié)果相符。

圖4 1980～2010年肥料施用總量和平均量

由圖5知，近些年來唐山、邯鄲、秦皇島等地化肥平均施用量提升較高，SOM含量增幅也較高；張家口耕地肥料施用提升最少，SOM含量增幅則較小。這表明耕地化肥施用量、耕作制度及種植結(jié)構(gòu)等社會經(jīng)濟活動與土壤肥力質(zhì)量變化總體上呈正相關(guān)，對SOM水平的提升具有一定正效應(yīng)。其中石家莊、保定、承德地區(qū)化肥施用增量居中，但其SOM增幅較小，可能是由于肥料利用率低、耕地管理差導(dǎo)致的，也可能是由于1980年時此處高值采樣點較多且過于集中，導(dǎo)致其代表性較差引起的。曾招兵等［27］分析得出廣東省隨著化肥的投入增加了作物產(chǎn)量，使得作物以根系或秸稈廢棄物形式歸還土壤的數(shù)量相應(yīng)增長，耕層土壤有機質(zhì)含量也相應(yīng)增加；趙吉霞等［28］研究表明1985～2006年云南墨江縣土壤全氮含量隨著氮肥施用量增加而產(chǎn)生了較小幅度的增加。上述對不同地區(qū)的研究結(jié)果表明，河北省不同地市SOM含量總體上可能隨化肥施用量增加而增加，但由于各地市消耗不同，具體增量不一樣，部分地市增量甚至為負(fù)值，導(dǎo)致河北省有機質(zhì)總量降低。

圖5 河北省各地市化肥平均施用增量與SOM平均增量關(guān)系

4 結(jié)論

從時間上看，1980～2010年間河北省SOM含量總體呈降低趨勢，平均降低3.6 g/kg。從空間上看，近30年全省SOM含量整體上均呈現(xiàn)西北減東南增的趨勢，增幅總體上由西北向東南依次增加；且其增加的面積占全省面積的65.54%，主要分布在南部平原地區(qū)和唐山、秦皇島等地。上述土壤有機質(zhì)變化格局的形成主要與施用化肥、實施秸稈還田等人類活動有關(guān)。因此，要想改善土壤肥力下降的問題，就應(yīng)堅持合理施用化肥、實施秸稈還田、制定合理的耕作制度、樹立用地養(yǎng)地相結(jié)合的耕作意識。

值得注意的是，本研究所采用的1980和2010年土壤樣點數(shù)據(jù)差距較大，而且空間分布位置和密度不一，所以通過回歸克里格插值模擬出的空間分布圖精度差距較大，其產(chǎn)生的誤差是否會對研究結(jié)果造成一定影響，也是值得探索的方面。