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

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

（華中農業大學資源與環境學院，農業農村部長江中下游耕地保育重點實驗室，湖北 武漢 430070）

土壤是人類進行各項生產和賴以生存的基地，土壤質量變化與質量高低關系到人類的生存狀況，而土壤有機質（SOM）是構成土壤有機無機復合膠體的核心物質，是土壤的重要組成部分，對土壤結構、物理性狀、化學性狀和生態過程有著重要的影響與作用［1-2］。所以，研究SOM的時空變化特征，對協調耕地土壤肥力，提高土壤質量水平，揭示土壤有機質空間分布格局，了解土壤質量動態變化，優化農業種植結構和施肥管理具有科學指導作用［3-4］。

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

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

1 材料與方法

1.1 研究區概況

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

1.2 數據來源

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

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

圖1 研究區域土壤樣點分布圖

1.3 研究方法

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

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

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

1.3.2 啞變量引入

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

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

1.4 分析工具

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

2 結果與分析

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

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

表1 第二次全國土壤普查土壤有機質分級標準 （g/kg）

表2 1980和2010年土壤樣點SOM描述性統計數據

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

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

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

2.2 回歸克里格插值結果

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

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

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

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

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

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

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

3 討論

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

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

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

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

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

4 結論

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

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