基于地統計學和空間自相關的土壤有機質空間異質性分析及分區*
——土地整治視角

左昕弘，賴佳鑫，劉 峰，石孝均

（西南大學資源環境學院，重慶 400715）

0 引言

耕地是土地資源的精華，是保障糧食安全的重要資源[1-2]。土地整治是保障我國糧食安全和生態安全的人類活動，也是耕地質量變化的重要驅動因素[3]。土地整治是指對低效利用、不合理利用、未利用及生產建設活動和自然災害損毀的土地進行整治、提高土地利用效率的活動[4-5]，農村土地整治更加強調通過田、水、路、林、村綜合整治，以增加有效耕地面積并提高耕地質量[6]。土地整治措施的實施對整治區域的土壤理化特性[7]、水資源空間分布[8]、景觀環境[9]、生物多樣性[10]等有重要影響[3]。如土地平整工程中田塊翻耕和表土剝離與回填等工程措施，通過增加土壤透氣性和保持土壤耕作層質量來改善土壤理化性質，從而提高耕地質量[11]。2021年4月通過的《中華人民共和國鄉村振興促進法》指出，縣級以上地方人民政府應當推進農村土地整治和農用地安全利用，加強農田水利基礎設施建設，改善農業生產條件。然而，未結合耕地土壤指標空間異質性的土地整治措施不僅不利于耕地質量改善，而且會增加整治成本[7]。因此，依據耕地質量指標空間異質性進行分區，并采取針對性的土地整治措施，是土地整治項目設計的重要基礎，對提高土地整治水平具有重要的意義。

土地整治分區是將實施土地整治工作的整個區域劃分為若干整治類型與方向相對一致的子區域。近年來，土地整治分區研究主要集中于研究尺度、評價指標和分區方法等方面。在研究尺度上，以宏觀區域（省[12]、市/縣[13]、自然區域[14]）居多，微觀（如項目區尺度[15]）的研究不足。在評價指標上，從自然質量、利用質量和經濟質量[16]逐步拓展到鄉村地域功能[17]、生態位[18]、生態系統服務[19]、生命共同體[5]、生態敏感性[20]等方面。但耕地質量指標，特別是土壤指標，仍是土地整治分區需要考慮的基礎和核心內容[21]。研究方法方面，在GIS空間分析的基礎上，多依據不同土地整治分區的側重點采用多因素綜合評價法[12]和屬性特征層次組合法[22]獲得分區，當研究對象事前分類方案不明確時，多采用聚類分析法[23]、空間自相關分析法[24]等進行分區。上述研究對提高土地整治效益，提升耕地質量提供了理論依據和技術支持，但因研究尺度較大，對土壤質量指標在土壤剖面上的異質性及水平方向上的空間自相關性少有考慮，從而忽略了評價指標的空間集聚程度及其顯著性，導致具有顯著空間集聚性的區域可能在分區時被割裂，不利于土地平整工程中表土剝離、移土培肥和土壤改良等關鍵技術的實施。

土壤有機質（soil organic matter，SOM）作為土壤質量和生產力最重要的指標之一，對改善土壤質地、提高土壤肥力、促進植物生長發育有重要作用[25]，也是基于耕地質量的土地整治分區中多選用的評價指標[26]。文章以土壤有機質為例，引入全局和局部空間自相關分析法，在渝西丘陵區典型土地整治項目區尺度上，探討基于土壤有機質垂直和水平空間上的異質性和空間自相關性的分區方法及其在土地整治分區中的應用。

1 研究區概況與研究方法

1.1 研究區概況

潼南國家農業科技園區位于渝西紫色土丘陵區。為建設耕地分布集中、基礎設施健全、農業發展條件較好的糧食主產區，從2017年底開始進行大規模的基本農田建設。該文研究區為2017年該園區基本農田集中建設區（105°44′39″E～105°45′26″E，30°04′03″N～30°04′53″N），在行政區劃上屬于重慶市潼南區柏梓鎮小嶺村，面積約120hm2（圖1）。研究區為淺丘帶壩地貌類型，海拔250～271 m；屬亞熱帶濕潤季風氣候，年均氣溫17.9℃，年均無霜期335d，年均降雨量975.5mm。土壤類型以侏羅紀中統遂寧組厚層泥巖母質發育的紅棕紫泥土和紅棕紫色水稻土為主，該土壤富含鈣質，呈微堿性。主要糧食作物為水稻和玉米。

1.2 土樣采集與測定

土壤樣品采集于土地整治前（2017 年5 月），依據研究區農用地地塊大小及形狀特征隨機布點采樣。研究區以河流（張家堰河）為界呈不同的地貌類型，河流以北為平壩區，農用地地塊面積較大，適當減少采樣密度，河流以南為丘陵區，受地形限制，農用地地塊較小，適當加大采樣密度。各采樣點按0～20cm、20～40cm、40～60 cm 進行分層采樣，分別表示為表層（Upper layer，U）、中層（Middle layer，M）和下層（Lower layer，L）。采樣時先去除地表凋落物，每個樣點向四周輻射共采集至少5 個分樣點，各土層土樣分別混合后四分法留取1 kg 土樣裝入樣品袋，共采集樣點53 個，分層土樣144 個（部分采樣點因土層薄，未采集到中、下層土樣）。同時用GPS記錄樣點經緯度和海拔（圖1）。

圖1 研究區地理位置、土地利用類型及采樣點分布

將采集的土壤樣品挑出礫石、草根等雜物，自然風干，研磨后在室內進行測定。SOM 含量采用重鉻酸鉀容量法—外加熱法測定[27]。

1.3 研究方法

1.3.1 地統計分析

普通Kriging 法（ordinary kriging）是地統計學中最常用的插值法，是建立在變異函數理論及結構分析基礎上，在有限區域內對區域化變量的取值進行無偏最優估計的一種方法。在GS+9.0中進行SOM 的半方差分析，以決定系數R2最大、殘差RSS最小確定最優理論插值模型。根據半方差擬合模型和參數，在Arc-GIS 10.2.2中采用普通Kriging法進行SOM的空間插值。

半方差函數的理論模型為[28]：

估算未采樣點的理論模型為[29-30]：

式（1）（2）中，?(h)為半方差函數，h為樣點空間距離，N(h)為空間距離為h的樣點數，r(xi)和r(xi+h)分別為變量r(xi)在空間位置xi和xi+h處的實測值，r(x0)為x0處的指標估計值，λi為第i個樣點的權重。

1.3.2 空間自相關分析

空間自相關分析（spatial autocorrelation analysis）用于檢驗具有空間位置的某要素的觀測值是否顯著地與其相鄰空間點上的觀測值相關聯[30-31]。包括全局空間自相關和局部空間自相關，通常采用Moran′s I指數表示空間要素的自相關程度。I 的取值范圍為-1～1，I＞0 表示變量在空間上呈正相關，表現出空間集聚特征，反之為負相關，表現出空間孤立特征，I=0 表示無空間自相關性。通常用Z值檢驗自相關顯著性，當Z＞1.96或Z＜-1.96時，變量的空間自相關顯著，否則不顯著，變量呈隨機分布。采用GeoDa 1.14進行Moran′s I分析，并在ArcGIS 10.2.2中結合Moran′s I分析結果制作LISA聚類圖。

全局Moran′s I指數的表達式為[30]：

局部Moran′s I指數的表達式為[30]：

式（4）中，N為空間變量個數，zi、zj分別為變量在i、j處的值，-z為變量z的均值Z值的計算公式為[32]：

其中：

式（6）中，Wij為zi和zj之間的空間權重函數，i≠j。

1.3.3 基礎數據處理

通過目視解譯及實地調查獲取2017 年研究區土地利用數據。運用Excel 2016 和SPSS 25.0對SOM 進行描述性統計分析和正態分布檢驗。采用ArcGIS10.2.2的Zonal Statistics工具生成各農用地斑塊的SOM含量。

2 結果與分析

2.1 SOM含量統計特征分析

研究區各土層SOM 含量的描述性統計特征如表1 所示。3 層SOM 的含量均值范圍為12.86～19.06 g/kg，表現為表層最高，且隨土層深度的增加而遞減。依據全國第二次土壤普查養分分級標準，各土層SOM 含量的均值均位于土壤養分分級標準的4 級。表層SOM 含量介于1～5 級，中、下層介于2～6 級（表2）。3 層SOM 的空間變異系數范圍為40.6%～60.8%，以表層最小，且隨土層深度的增加而逐漸增加，均屬于中等變異。

表1 不同土層SOM的描述性統計特征

表2 中國土壤養分分級標準（全國第二次土壤普查）

2.2 SOM的空間變異性分析

經典統計學方法僅在一定程度上反映研究區SOM 的樣本總體特征及變異狀況，不能定量地刻畫其隨機性和結構性、獨立性和相關性，需進一步采用地統計學方法對其空間變異結構進行分析和探討[33]。經K-S 檢驗，3 層SOM 除表層呈正態分布外，其余土層經對數轉換后呈正態分布，滿足地統計學分析要求。該文采用GS+9.0對不同土層SOM 的空間變異性進行分析，用球狀、指數和高斯等模型進行擬合，得到模型的相關參數值，當決定系數（R2）接近于1且殘差（RSS）較小時，模型擬合度最優。如表3所示，3層土壤中，除中層SOM 的最優擬合模型為球面模型外，其余均為高斯模型。最優擬合模型結構參數中塊金值（C0）反映的是最小抽樣尺度下變量的測量誤差和隨機因素引起的變異。基臺值（C0+C）表示系統內的總變異量，基臺值較大表明系統變量具有較強的空間變異。塊基比C0/（C0+C）表示系統變量的空間相關性程度，若比值＞25%，表明系統具有強烈的空間相關性，若比值介于25%～75%，表明系統具有中等的空間相關性，若比值＞75%，表明系統空間相關性較弱，若比值接近于1，表明在整個尺度上具有恒定的變異[34-35]。研究區3 層SOM 的C0/（C0+C）介于28.6%～33.3%，具有中等的空間相關性，表明各土層SOM 的空間變異受結構性因素（成土母質、土壤類型、地形地貌等）和隨機性因素（耕作、施肥管理、種植制度等）的共同影響。變程指空間相關性的作用范圍，反映區域化變量影響范圍的大小，同時也反映變量自相關范圍的大小。研究區SOM的變程介于88.33～102.19 m，表明SOM含量在該范圍內存在空間相關性。

表3 不同土層SOM的半方差函數理論模型和相關參數

2.3 SOM的空間分布特征分析

為更直觀地揭示各土層SOM 在研究區內的空間分布格局，該文基于最優半方差擬合模型所得參數，采用普通Kriging 法在ArcGIS 中進行插值，得到3 層SOM 的空間分布圖。如圖2 所示，各土層SOM 呈不規則的斑狀與塊狀分布。SOM 含量在表層表現為西北部和東南部的部分區域較高，中部和北部較低；在中、下層表現為由北向南遞減。此外，表層和中層SOM 含量介于10～20 g/kg之間的區域面積最大，下層則以＞10 g/kg的區域為主。

2.4 SOM的空間自相關分析

2.4.1 全局空間自相關分析

為計算研究區農用地斑塊的全局與局部Moran′s I指數并繪制散點圖及LISA聚類圖，該文基于SOM空間分布圖（圖2），采用ArcGIS的Zonal Statistics工具，生成各農用地斑塊的SOM平均含量作為基礎圖層（圖3）。

圖2 不同土層SOM的空間分布

圖3 不同土層農用地斑塊SOM平均含量的空間分布

為描述SOM 的空間自相關性隨滯后距的變化規律，該文采用GeoDa 軟件和Excel 計算并繪制了各土層SOM 在不同滯后距下的空間自相關程度。如圖4 所示，各土層SOM 的全局Moran′s I 指數在滯后距位于68.03～1 000.00 m 的范圍內大于0，表現為較強的正相關性，表明SOM 在上述范圍內呈正空間自相關，存在明顯的空間集聚性。當滯后距為68.03 m 時，在保證每個對象都有至少1 個鄰居的情況下，各土層SOM 的Moran′s I 指數均達到峰值，此時SOM 的正空間自相關集聚性最強。隨著滯后距的增加，Moran′s I 指數下降，且下降趨勢逐漸減緩，表明以68.03 m 為閾值構建的空間權重和計算得到的空間自相關結果具有更高的精確度。因此該文選擇68.03 m 為閾值對SOM進行空間自相關分析。

圖4 不同土層SOM的空間自相關性隨滯后距的變化

為探討各土層SOM 的空間關聯程度，該文以農用地斑塊為基本單元，對各土層SOM 進行全局空間自相關分析，計算全局Moran′s I 指數并繪制Moran′s I 散點圖。如表4 所示，各土層SOM 的Moran′s I指數均大于0.8，呈較強的正空間自相關集聚性，表明各土層SOM 的空間分布表現為相似值的空間集聚，即高值趨向于與高值相鄰，低值趨向于與低值相鄰，且相鄰觀測值間相似性較強，空間變異性則相對較弱。各土層SOM的Z值均大于1.96，呈顯著空間自相關，表明觀測值與周圍觀測值相似。顯著性檢驗P值均小于0.01，表明各土層SOM含量均表現為較強的正空間自相關集聚態勢。

表4 不同土層SOM的空間自相關參數

Moran′s I散點圖中第一、三象限為高—高（HH）和低—低（LL）集聚區域，表明指標存在正的空間關聯，HH 和LL 分別表示高值斑塊被其他高值斑塊所包圍和低值斑塊被其他低值斑塊所包圍；第二、四象限為低—高（LH）和高—低（HL）集聚區域，表明指標存在負的空間關聯，LH 和HL 分別表示低值斑塊被其他高值斑塊所包圍和高值斑塊被其他低值斑塊所包圍。如圖5 所示，各土層的SOM 含量主要分布于一、三象限，表現為明顯的HH或LL集聚特征，表明研究區SOM的空間自相關性整體呈現為兩大趨勢，即SOM 含量高值區較集中，形成高值與高值自相關趨勢，同時，低值區也較為集中，形成低值與低值自相關趨勢。

圖5 不同土層SOM的全局Moran's I散點

2.4.2 局部空間自相關分析

全局空間自相關的Moran′s I 指數可判斷各指標是否存在空間集聚區及空間離散區，但無法獲得其在空間上的分布狀況，而局部空間自相關指標結合LISA 聚類圖能夠揭示各指標在空間上具體的分布位置。因此，為進一步明確各土層SOM 含量空間集聚區和空間離散區的區域分布特征，采用局部空間自相關指標及LISA 聚類圖將農用地斑塊劃分為HH、HL、LH、LL和“不顯著”（NN）5種類型并繪制空間分布圖。各類型斑塊數量與面積分別占總斑塊數量與面積的比例如表5所示，SOMM表現為HH型占顯著性斑塊的比例最大，表明中層SOM主要表現為高—高集聚特征；其余均為LL型所占比例最大，即表、下層SOM主要表現為低—低集聚特征，也表明在土地整治中保護表層養分含量較高地塊的重要性。此外，各土層SOM均無LH型，表明研究區無SOM低值被高值所包圍的農用地斑塊。

表5 各類型斑塊數量與面積占比 %

LISA 聚類圖表明（圖6），表層SOM 的HH 集聚區主要分布于平壩區西北部和東南部丘間谷地，LL 集聚區主要分布于南部丘陵區，部分位于平壩區西部；中、下層SOM 的HH 集聚區則主要集中于北部平壩區，LL集聚區主要分布于南部丘陵區。

圖6 不同土層SOM的LISA聚類

2.5 SOM分區

基于SOM 空間自相關分析結果，依據以下原則對耕地進行分區：（1）以表層SOM 空間自相關分布特征為主，中、下層為輔，將3 層SOM 的高值集聚分布、低值集聚分布、高低離散分布、低高離散分布和空間隨機分布疊加組合進行分區；（2）3層SOM空間自相關類型疊加后，將斑塊面積較小的類別根據表層SOM 空間自相關類型合并到其他斑塊面積較大的分區。結果表明，研究區可劃分為高SOM 集聚區、表層高SOM集聚區、低SOM集聚區和SOM隨機分布區4個分區（表6），空間分布如圖7所示。

高SOM集聚區主要分布在水稻種植區，共有31個耕地斑塊，面積為4.98hm2，占研究區耕地總面積的5.2%，該區域宜作為優先整治區域（表6 和圖7）。由圖8 可知，該區域表層SOM 含量顯著高于中、下層（P＞0.05），3 層SOM 含量均顯著高于其他分區（P＞0.05），是各土層SOM 含量最高的耕地集聚區。該區域自然條件較好，耕地集中連片分布，為耕地利用提供了良好的基礎。在土地整治中應重視表土剝離與回填，同時應保護和維持現有自然條件和利用管理水平，避免經濟建設對耕地質量的影響。

表層高SOM 集聚區主要分布在東南部丘間谷地，共有31 個耕地斑塊，面積為5.43 hm2，占研究區耕地總面積的5.6%，該區域宜作為適度整治區（表6 和圖7）。由圖8 可知，該區域表層SOM 含量顯著高于中、下層（P＞0.05），且表層SOM 含量顯著高于除高SOM 集聚區外的其他分區（P＞0.05），中層SOM 含量顯著高于低SOM 集聚區（P＞0.05），下層SOM 含量略高于低SOM 集聚區，是表層SOM 含量相對較高的耕地集聚區。在土地整治中應同樣重視對表層的保護，同時應加強對中、下層SOM 含量的提升（如有機物料回填、機械改土等）。此外，應保護和維持現有自然條件和利用管理水平，盡量避免經濟建設對耕地質量的影響。

圖7 SOM分區

低SOM 集聚區主要分布在平壩西部和南部丘陵區坡耕地區域，共有111 個斑塊，面積為14.56hm2，占研究區耕地總面積的15.1%，該區域宜作為重點整治區（表6 和圖7）。由圖8 可知，該區域表層SOM 含量顯著高于中、下層（P＞0.05），3層SOM 含量均低于其他分區，且除下層外均表現出顯著性（P＞0.05），是各層SOM 含量最低的耕地集聚區。在土地整治中應結合自然條件和利用管理特征進行中長期改良措施（如坡改梯、改良培肥等），同時改善農田生產條件，以提升區域整體耕地質量水平。

SOM 隨機分布區空間分布較為廣泛，共有327 個斑塊，面積為71.53 hm2，占研究區耕地總面積的74.1%，該區域宜作為優化調整區（表6 和圖7）。由圖8 可知，該區域表層SOM 含量顯著高于中、下層（P＞0.05），3 層SOM 含量均處于中等水平。在土地整治中需注意對少量優質耕地的保護，以此為中心優化農田布局，同時應考慮適當調整、優化農業結構。

圖8 各分區的SOM含量大寫字母表示同一分區下SOM在不同土層間的顯著性差異，小寫字母表示同一土層下SOM在不同分區間的顯著性差異；誤差線代表±SD

表6 SOM分區及區域特征

3 討論

半方差函數分析結果表明，研究區各土層SOM 具有中等的空間相關性，受隨機性因素與結構性因素的共同影響，空間自相關范圍較小。滯后距較小時，SOM 表現為較強的正空間自相關性，存在明顯的空間集聚性（圖4 和表4）。這與半方差函數分析的結果較為一致（表3）。但二者對空間自相關衡量的嚴格程度不同，半方差函數的分析結果包括了正負自相關，而全局Moran′s I 指數將正負自相關進行了區分[36]。從空間分布上看，各土層SOM 呈不規則的斑狀與塊狀分布（圖2），LISA 聚類圖也表現為圍繞HH 和LL的多中心分布格局（圖6），這一結果可以指導SOM 在分區時的精準劃分，進一步為土地整治分區提供依據。

結合土地利用類型來看，研究區水田所處地勢平緩，積水條件下SOM 分解較緩慢，易于積累，且水田的秸稈還田率較高，被微生物轉化為有機碳的量較多[37]，使得水田表層SOM 最高且顯著高于撂荒地和旱地（P＞0.05）（圖9），因此SOM 表層的HH 集聚區主要為水田（圖1），分布在北部平壩區和東南部丘間谷地。受施肥和表層土壤腐殖質隨下滲水在土體中淋溶、遷移、淀積的影響[38]，中、下層SOM 含量也以水田最高且顯著高于其他土地利用類型（P＞0.05），其次為旱地，因此中、下層SOM 的HH 集聚區除水田外還有部分旱地，主要分布在北部平壩區。這是土地利用類型影響土壤養分空間變異及剖面異質性的直接體現。但LL 集聚區土地利用類型表現多樣化，且主要分布在南部丘陵區，表明除土地利用類型外，SOM 還受到地形因子、成土母質、土壤類型和秸稈還田方式等其他因素的影響[39]。

圖9 不同土地利用類型下不同土層的SOM含量大寫字母表示同一土地利用類，SOM在不同土層間的顯著性差異，小寫字母表示同一土層下SOM在不同土地利用類型間的顯著性差異；誤差線代表±SD

基于SOM 的剖面異質性和空間自相關性劃分的分區，可服務于土地整治中的均質單元劃分、表土剝離與回填、移土培肥和土壤改良等工程措施的實施。其中，全局空間自相關明確了SOM 的空間相關性與規律性，局部空間自相關進一步描述了SOM 的空間自相關程度，Moran′s I 散點圖和LISA 聚類圖則展示了上述空間自相關在空間上的具體分布格局[40]。與依據SOM 含量等級劃分的分區相比，基于3 層SOM 的LISA 聚類結果的疊加分區，能將SOM 在垂直和水平方向上的分布情況與空間集聚性結合，從而提高評價指標分區及基于此的土地整治分區的合理性，更有利于針對性地采取土地整治措施。

4 結論

在渝西典型紫色土丘陵區，0～60 cm 土層中的SOM 含量總體較低，受結構性和隨機性因素的共同影響，各土層SOM 均表現為顯著的正空間自相關，且空間自相關范圍較小。在水平空間分布上，各土層SOM 呈不規則的斑狀與塊狀分布，LISA 聚類圖也表現為圍繞HH 和LL 的多中心分布格局。土地利用類型對SOM 的空間分布存在顯著影響，HH 集聚區主要為水田，多分布于平壩和東南部丘間谷地，LL 集聚區土地利用類型則表現多樣化，主要分布于南部丘陵區。以表層SOM 空間自相關分布特征為主，中、下層為輔，可將研究區劃分為高SOM 集聚區、表層高SOM 集聚區、低SOM 集聚區和SOM 隨機分布區4 個分區。各分區體現了SOM 含量在剖面和水平上高值與低值的集聚性特征，可作為土地整治分區，尤其是土地平整工程中的表土剝離、移土培肥和土壤改良等關鍵技術實施的重要參考。