平原丘陵過渡區土壤有機質空間變異及其影響因素

杜佩穎 張海濤 郭 龍 楊順華 章 清 田 雪

（華中農業大學資源與環境學院，武漢 430070）

土壤是在成土因素的長期相互作用下形成的復雜自然綜合體，具有高度的空間異質性和依賴性。土壤有機質（Soil Organic Matter, SOM）是土壤肥力的主要來源，對農業生產和可持續發展具有重要意義，因此了解土壤有機質的空間分布，掌握其空間變異規律對生產實踐具有重要的指導意義［1-3］。

目前地統計模型已經被成熟地應用于土壤有機質空間變異規律的研究中，以變異函數為工具的區域化變量理論得到了廣泛應用［4-5］。普通克里格方法（Ordinary Kriging, OK）作為經典的地統計模型較好地考慮了土壤樣點屬性的空間異質性和依賴性，但OK未將環境影響因子納入模型，忽視了土壤屬性在形成、發育和遷移過程中會受到的多種自然社會因素的影響。因此在構建土壤屬性模型時考慮不同環境因素對土壤有機質空間變異規律的影響具有重要意義，為建立精確的土壤屬性空間模型提供有價值的輔助信息［6-7］。

眾多研究表明土壤有機質受到氣候、植被、地形、土壤營養元素、土壤質地、土壤pH以及人為活動等因素的綜合影響［8-10］。因此為進一步準確描述土壤有機質及其影響因素之間的響應機制，學者們將土壤有機質的異質性及其環境影響因素納入模型，進行了土壤有機質預測和制圖等方面的研究：例如，連綱等［7］以黃土丘陵溝壑區縣域土壤有機質為研究對象，利用多元線性回歸模型進行空間預測，結果表明包含地形濕度指數、高程、林地、梯田和土壤調節植被指數的區域回歸模型具有最優的預測效果；郭月峰等［11］利用老哈河流域的地形因子對土壤有機碳建立多元線性回歸模型，結果表明高程和坡度對該區域土壤有機碳影響作用較強。這些研究結果不僅進一步表明了環境變量與土壤有機質之間的密切聯系，更為進一步進行土壤有機質預測及制圖提供了發展方向及理論基礎。

多元線性回歸模型將整個研究區域當作一個研究整體，未考慮回歸系數的空間局部差異，因此線性回歸模型雖然可以揭示土壤有機質與環境變量之間的關系，但是難以建立土壤有機質與環境變量關系在局部區域上的差異性［12-13］。地理加權回歸（Geography Weighting Regression，GWR）模型作為局部線性回歸模型，不僅考慮了土壤屬性的空間異質性，同時考慮了環境變量與土壤屬性之間的相關關系的空間非平穩性［14-16］。GWR模型可以根據土壤屬性的空間變異規律調節并選擇合適的帶寬，反映環境影響因子對局部地理位置土壤有機質的影響程度。GWR模型目前較多地被應用于土壤屬性的預測和模擬中，但是由于GWR模型的前提假設是納入模型的變量和殘差具有獨立性，因此土壤有機質和環境變量在空間上的依賴性會違背GWR模型的前提假設，并且影響模型性能的可靠性，為此GWRK模型由于其考慮了回歸殘差自相關性的獨特優勢被廣泛地應用于土壤學領域［24, 26］。

本研究區內土壤有機質影響因子復雜，土壤有機質在局部區域表現出很強的異質性，因此利用傳統的地統計模型難以準確模擬土壤有機質的空間局部特征，而多元線性回歸模型無法考慮環境變量與土壤有機質的空間異質性［21］。為此本文選取GWR模型作為基礎研究模型，探討影響因子對土壤有機質局部區域影響作用的差異性，同時利用GWRK模型將GWR模型的預測結果與OK進一步結合，以達到更好的預測精度和效果。楊順華等［20］曾以此區域為研究對象，將地形因子納入回歸模型對土壤有機質含量進行了預測，并得到GWRK具有較高預測精度的結論，但是該研究沒有進一步探討土壤有機質與自身物理化學性質之間的關系，為此本文進一步完善輔助數據集，為過渡地帶復雜環境下數字土壤制圖提供研究基礎，同時該研究也未進行土壤有機質影響因子空間非平穩性影響的研究，本文通過影響因子系數分布圖深入分析了各影響因子在局部區域對土壤有機質影響作用的變化規律。本文的主要研究內容為（1）建立由地形因子和土壤屬性因子組成的土壤有機質預測輔助數據集，（2）運用GWRK預測土壤有機質含量并得到其空間變異規律，為平原丘陵過渡地帶數字土壤制圖提供參考依據，（3）分析土壤有機質影響因子變化規律，結合研究區地形特點和土壤有機質分布情況分析影響因子對土壤有機質影響作用的空間非平穩性。

1 材料與方法

1.1 研究區概況

本研究區位于宜都市枝城鎮境內（圖1）。宜都市位于3 0°0 5′～3 0°3 6′N，111°05′～111°36′E之間，地處湖北省西南部，江漢平原西部。枝城鎮位于宜都市東北部，是江漢平原向鄂西山區的過渡地帶，亞熱帶季風氣候，四季分明，雨熱同期，年平均降水量1 212 mm，氣溫16.7℃，無霜期273 d，地形地貌以平原、丘陵為主，地勢西南高，東北低，大部地區屬武陵山脈的丘陵地帶，東北部臨長江有部分沖積平原。結合地形地貌，在實地調查的基礎上選取一塊面積100 km2地形起伏較明顯的代表性區域作為研究區，研究區最高點海拔為483 m，最低點為0 m，高程起伏變化較大，土壤類型主要有黃棕壤土類、石灰巖性土和水稻土類等，土地利用類型主要有林地、水田、旱地和園地，水田的種植方式為水稻輪作，一年兩季，旱地多種植玉米、土豆等農作物，園地主要種植柑橘和茶樹。

1.2 樣品采集與分析

500個樣點土壤樣品采集于2013年12月—2014年1月，布點方法采用“網格+隨機”相結合，遵循平原區域較稀、山區加密的原則，采樣方法主要使用“梅花法”，對于較狹窄的地帶或農田采用“蛇形法”。土樣選取0～20 cm表層土壤，并使用差分式全球定位系統（DGPS）記錄實地樣點的空間位置。土壤有機質含量采用重鉻酸鉀氧化—外加熱法測定，機械組成采用馬爾文—激光粒度儀測定，使用美國制分級：砂粒（2～0.05 mm），粉粒（0.05～0.002 mm），黏粒（＜0.002 mm），土壤容重采用環刀法測定［22］。

圖1 研究區位置及采樣點分布（392個建模點和98個驗證點）Fig. 1 Location of the study area and distribution of sampling points (including 392 for modeling and 98 for validation in the study area

1.3 輔助信息獲取

研究區的數字高程模型（D E M）數據（30 m×30 m）來源于國際科學數據共享平臺（http:/ /datamirror.csdb.cn /admin /datademMain.jsp），利用DEM數據獲取其衍生數據變量，主要為坡度（Slope）、坡向（Aspect）、高程（Elevation）、粗糙度（Roughness）、地形濕度指數（Topographic Wetness Index，TWI）、匯流動力指數（Stream Power Index，SPI）、沉積物運移指數（Sediment Transport Index，STI）。其中坡度采用ArcGIS10.2中空間分析模塊中的相關工具直接提取，其余環境因子使用ArcGIS10.2中的柵格計算器以及水文分析模塊綜合計算獲取。

1.4 提取影響土壤有機質含量的重要指標

除地形因子外，土壤屬性因子也是影響土壤有機質含量的重要指標，因此除上文提到的地形因子外，加入有效鐵(AI)、土壤pH、容重（BD）、礫石度（GD）、歸一化植被指數（NDVI）、亞鐵礦物指數（FMI）和土壤機械組成（包括黏粒、粉粒、砂粒）這7個屬性作為輔助變量，組成由14個輔助因子組成的初選因子集，將土壤有機質含量與影響因子進行皮爾遜相關分析和逐步回歸分析，篩選出與土壤有機質含量顯著相關的影響因子。利用逐步回歸方法選取進入模型的回歸變量，既能保證與土壤有機質含量極顯著相關的輔助因子進入回歸模型，又能有效去除自變量之間的共線性，使用該方法得到研究區GWR模型的最佳解釋變量：Elevation、Slope、Aspect、AI、BD、GD、Clay。

1.5 模型方法

（1）回歸克里格（RK）方法［25-26］首先對解釋變量和目標變量之間的回歸關系進行探討，建立兩者之間的線性回歸關系，進而由回歸關系得到目標點處的確定性趨勢項，對分離趨勢項后的殘差進行普通克里格（OK）插值得到殘差項，最后將趨勢項和殘差項相加得到采樣點處的RK預測值。

（2）地理加權回歸克里格（GWRK）［12,19,27］是一種局部回歸方法，它將全局回歸方法中的全局回歸值換成GWR局部回歸值，同時考慮了局部預測的殘差。GWRK是對GWR的延伸與擴展，即對局部模型GWR擬合后得到的殘差進行OK插值，然后與GWR擬合的趨勢相加。

（3）數據標準化。在評估解釋變量對SOM的影響時，預測模型中的系數是一個可以定量比較的指標，但在進行這一過程之前，必須對解釋變量進行標準化轉換，將解釋變量的量綱進行統一，標準化后的變量其系數越大說明它對土壤有機質的影響程度越高，標準化公式如下：

式中，χ*為標準化結果，χ 為解釋變量的原始值，μ為變量的均值，σ為變量的標準差。

（4）模型精度驗證。為綜合比較模型的預測精度，從數值上精確比較預測效果，本文選擇平均誤差（ME）、平均絕對誤差（MAE）、均方根誤差（RMSE）、皮爾遜相關系數（Pearson’r）以及不精確度（IP）等指標評估預測精度。本研究利用箱線圖法剔除異常值10個，利用ArcGIS中的工具將490個樣點隨機分為兩組，其中一組392個（占樣點總數的80%）作為建立模型的擬合數據集，另一組為其余98個（占樣點總數的20%）作為驗證數據集用于驗證模型。

2 結 果

2.1 土壤有機質描述性統計特征

由表1可知，建模點的土壤有機質含量在3.80～69.40 g·kg-1之間，平均值為28.42 g·kg-1，變異系數為39.86%，屬中等變異，土壤有機質含量不符合正態分布，經對數轉換后通過K-S檢驗，符合正態分布。研究區高程均值為172.3 m，最小值為0 m，最大值為483 m，變異系數為65.06%，具有較強的變異程度，高程變化幅度較大；坡度和坡向的均值分別為9.31°和177.2，變異系數分別為80.02%、58.49%，變異系數較大，說明地形起伏程度較大；有效鐵、土壤礫石度、黏粒的變異系數分別為93.89%、162.1%和51.10%，研究區土壤屬性因子變異程度較大。

表1 土壤有機質含量及其影響因子的描述性統計結果Table 1 Descriptive statistic of soil organic matter (SOM) and its affecting factors

2.2 土壤有機質預測模型的建立

土壤有機質預測模型中使用高程、坡度、坡向、有效鐵、容重、土壤礫石度和黏粒7個影響因素作為解釋變量建立多元線性回歸模型（MLR）和地理加權回歸（GWR）。MLR模型參數的統計結果表明，高程、坡度和土壤容重系數的絕對值大于坡向、有效鐵、土壤礫石度和黏粒含量，系數越大表明影響程度越大。土壤有機質含量與高程、坡度、坡向、有效鐵和土壤礫石度呈正相關，而土壤有機質與土壤容重和黏粒呈負相關。方差膨脹因子（VIF）在逐步線性回歸中反映了冗余的解釋變量信息，VIF大于7.5的解釋變量存在局部共線性需要去除，本模型中VIF均小于7.5，7個解釋變量均適合納入模型。

GWR是局部回歸模型，解釋變量的回歸系數在研究區內隨空間位置的不同隨之變化， GWR模型參數的統計結果顯示，解釋變量的范圍、最小值、最大值、平均值、標準差和變異系數如表所示。GWR模型中的參數與MLR有著類似的含義，系數范圍、標準差和變異系數用來描述解釋變量空間系數的分布情況，除土壤容重外，其他解釋變量均具有強空間變異性（＞35%）。

2.3 土壤有機質預測模型

由表2可知，土壤有機質含量、線性回歸殘差和地理加權回歸殘差的理論半方差模型分別為指數模型、指數模型和球狀模型。塊金值分別為0.021、0.032、0.087，說明土壤有機質存在較小的隨機性誤差。塊基比分別為0.23、0.19、0.16，說明三個模型存在不同程度的空間自相關性，且均小于25%，適合進行普通克里格插值。變程分別為142 3 m、134 3 m和158 0 m，在此范圍內的空間變量具有空間自相關性，變程以外則不存在；三者的決定系數均在80%以上，取得了較好的模型擬合效果，較小的殘差平方和也表明了擬合點和擬合曲線之間較好的吻合效果。

表2 有機質、線性回歸殘差和GWR殘差的半方差函數模型及參數Table 2 Semi-variance models of SOM, OLS and the residuals of GWR and parameters involved

采用OK、RK和GWRK三種方法分別預測研究區土壤有機質含量（圖2），從圖中可以看出，GWRK預測得到的SOM區間范圍較OK和RK更大；三種方法的預測結果在整體上趨勢一致，高值區集中分布在研究區東北部，低值區聚集在南偏東方向，在研究區南部及西南部也有部分低值區域存在；GWRK在低值區域的預測范圍更大，OK和RK在高值區域的預測范圍更大，整體上呈現出“高低值區域縮小，中值區域擴大”的趨勢；在高低值過渡的區域，OK和RK的過渡范圍小，過渡幅度較大，存在過渡區突變的情況，過渡不夠平緩，而GWRK在過渡方面則體現出更好的效果，過渡較為平緩，整體上沒有較突兀的顯示效果，高低值界線更加模糊化，過渡曲線更加曲折，這與實際SOM分布規律更加吻合。從空間分布圖上來看，GWRK具有更好的預測效果，這主要是因為GWRK更好地考慮了局部影響因子對土壤有機質的影響作用，更加全面、準確地將解釋變量納入模型，同時考慮了預測殘差，有利于反映研究區細節變化特征，使預測結果與真實值更加接近，有助于掌握SOM的整體分布情況。

2.4 預測模型精度評價

為進一步說明GWRK預測精度，引入模型評價指標對預測模型進行精度評價。選用平均誤差（ME）、平均絕對誤差(MAE)、均方根誤差(RMSE)、不精確度（IP）和相關系數（r）五個指標進行模型評價。從表3可以看出，三者的ME均為正值，說明預測值整體水平較實測值偏低，GWRK平均誤差最小；GWRK的MAE和RMSE均小于其他兩者的結果，降低程度分別為26.76%和39.28%；不精確度則是OK、RK高于GWRK；GWRK的相關系數 r 高于OK和RK。由定量數據分析可知，局部預測模型GWRK較OK和全局模型RK預測精度更高。

2.5 影響因子的空間非平穩性分析

為說明解釋變量對SOM影響作用的變化規律，繪制解釋變量系數空間分布圖，如圖3。

圖2 土壤有機質含量空間分布圖（a.普通克里格，b.回歸克里格，c.地理加權回歸克里格）Fig. 2 SOM spatial distribution map interpolated with OK, RK and GWRK

表3 普通克里格、回歸克里格和地理加權回歸克里格的精度評價Table 3 Comparison between OK, RK and GWRK in precision

研究區地勢西高東低，為典型的平原丘陵過渡地帶，從高程系數分布圖(圖3b)可以看出，高程對東北部影響程度最弱，這可能與東北部較為平坦的地勢有關；從東北向西南遞進，出現了高程對SOM影響程度最大的區域，因為此處正是研究區平原丘陵地貌分界線（200 m）所在的區域，因此也是影響程度最高的區域。坡度的不同引起地表徑流強弱的不同，同時坡度較大的地方土地利用類型為林地的可能性最大，地表大量枯枝落葉阻礙了水分的下滲淋溶，樹木對徑流的緩沖作用也有利于SOM積累［18］，由坡度系數分布圖(圖3c)可知，西部和東部受坡度影響作用最大，從兩側向中間移動，影響作用逐漸減小，最低值區分布在研究區北部邊緣，坡度越小，SOM與坡度空間變化一致性越強，反之隨坡度增加，SOM影響因素復雜化，因此坡度對SOM的空間影響權重減小。坡向決定了區域接受太陽輻射的角度和強度，影響光熱資源再分配，接受太陽輻射多的區域有機質循環速度快，更有利于SOM積累，研究區坡向的正弦值和余弦值分別代表坡向朝東和朝北的程度，且前者與SOM呈負相關，后者呈正相關，表明研究區西北坡向更有利于有機質的積累，坡向系數分布圖(圖3d)顯示，研究區西北部受坡向影響最大，且從西北部向中部影響逐漸減弱，這與太陽輻射規律吻合。有效鐵系數分布圖(圖3e)(AI)顯示，研究區中西部受有效鐵影響作用最大，這可能與過渡區高程變化有關，此處向南北兩側影響作用減弱，至東南角處減為最弱，而東北部地勢較平坦處受有效鐵影響作用中等，且強弱分布較均勻，這是因為在地勢平坦處有效鐵隨徑流變化較小，含量分布均勻。容重系數分布圖(圖3f)和礫石度系數分布圖(圖3g)的強弱變化規律相似，趨勢相反，說明兩者對SOM的影響作用相反，容重系數從中心至四周逐漸增大，最小值出現在中心位置，最大值位于研究區邊緣，而礫石度系數則是從中心至四周逐漸減小，分布規律相反，由此可知，受土壤容重影響作用強的地方土壤礫石度影響作用弱，呈負相關關系。土壤機械組成［13］（黏粒、粉粒、砂粒）通過土壤的通透性、保水保肥性等影響SOM，尤其是砂粒和黏粒在局部范圍對SOM影響作用較為顯著，黏粒會影響土壤物理和化學保護機制，黏粒含量越高，土壤對SOM的吸附作用更強，降低SOM礦化速度，有利于SOM積累；砂粒的影響機制恰恰相反，砂粒含量較高的土壤，有機質礦化速度加快，不利于有機質的積累，本研究中黏粒（Clay）對SOM影響作用顯著，從黏粒系數分布圖(圖3h)可以看出，研究區自西北至東南影響作用逐漸減弱，最高值區位于中西部區域，最低值區位于東南部，東北地勢平坦區域影響作用變化微弱。

圖3 解釋變量系數的空間分布圖Fig. 3 Spatial distribution of explanatory variable coefficients

3 結 論

土壤有機質是影響土壤肥力的重要指標，掌握其影響因素的變化規律對指導農業生產具有重要意義。本文所選研究區為平原丘陵過渡地帶，地形條件復雜，人為活動頻繁，運用GWRK對建模集392個采樣點建模預測，通過分析得到研究區模型最佳解釋變量為：高程、坡度、坡向、有效鐵、容重、礫石度、黏粒，且GWRK由于考慮了區域化變量的局部非平穩性，較OK和RK具有更好的預測效果。計算各影響因素權重，并繪制其系數分布圖，從圖中可以得到各影響因子在不同空間位置對土壤有機質含量的影響和變化規律，從而針對性地采取農業改造措施，在影響因素作用較強的區域重點開展農業設施或局地改造，充分利用地形特點和土壤屬性特征，合理開展農業布局，因地制宜地進行農業生產，為實踐生產提供參考依據。