基于GIS的環巢湖地區土壤養分空間變異特征研究

童童，梅帥，劉瑩，常珺楓，梁華忠，范其龍，王強，馬友華*

（1.安徽農業大學資源與環境學院，合肥 230036；2.合肥市農業經濟技術服務管理總站，合肥 230000）

糧食是人類賴以生存的資本，隨著人類的發展，糧食的需求量也在不斷提高。施肥是保障糧食產量的重要手段，但盲目施肥導致養分流失嚴重[1]，造成河流湖泊等水體污染問題。根據土壤養分狀況，合理配施有機無機肥、選擇適宜的種植模式，能夠有效提高肥料利用率，改善土壤理化性質，減少農業面源污染。土壤養分的空間變異由母質、地形、氣象等結構性因素和施肥、種植習慣等人為活動的隨機性因素共同影響[2]，研究土壤養分的空間變異特征對充分掌握養分分布規律及促進農業可持續發展具有重要的現實意義。

地統計學中的半變異函數法和Moran's I 指數法已被廣泛地應用于空間自相關性的研究[3]。已有學者[4]通過在研究區域選取土壤樣點并檢測其養分指標，利用地統計學方法和GIS 技術探討了土壤養分在研究區域的分布特征和空間異質性；陳文舉等[5]則通過克里金插值法和單因素方差分析，研究了不同深度土壤pH 的空間分布特征，并了解了土壤酸化分層現象及酸化趨勢；Chen 等[6]獲取了不同年份的土壤樣點數據，通過半方差函數法和克里金插值法，研究了土壤養分的時空變化特征。半變異函數的參數能夠表明其空間自相關程度，而Moran's I 指數能夠反映出研究對象的分布聚集狀況。王強等[7]利用遙感影像數據提取出不同種植年限和種植模式下菜地的空間分布，并以1 km 的網格為評價單元計算局部Moran's I 指數，從而識別出菜地土壤養分空間聚類和空間異常值，最終結果表明菜地土壤養分受到種植時間與城鎮距離因素的影響；陳清霞等[8]將半變異函數與Moran's I指數結合，分析出土壤pH 具有小范圍、中等程度的空間自相關性，且其空間自相關性顯著，體現在土壤pH呈斑塊狀的鑲嵌分布。大多數針對土壤養分空間變異特征的研究是在小范圍尺度或大范圍尺度低密度采樣的數據前提下，其研究結果的精度受限，且研究的土壤養分指標較為單一，少有養分之間全面的對比研究。

本研究以巢湖周邊地區（巢湖市、肥東縣、肥西縣、廬江縣、包河區）為研究范圍，以同一時期高密度采樣下土壤養分（pH、有機質、全氮、堿解氮、有效磷、速效鉀、緩效鉀）檢測結果作為研究數據。基于ArcMap10.6軟件，計算并選擇最佳半變異函數參數擬合模型，分析其空間自相關程度；結合全局Moran's I指數反映土壤養分自相關性的大小及其是否顯著[9]；在最佳擬合模型的基礎上，用普通克里金插值法得到養分分布圖，并計算局部Moran's I 指數制作LISA 聚類圖。通過將半變異函數與全局/局部Moran's I 指數結合，從而得到更加全面有效的研究結果[10]。本研究結果可為環巢湖地區測土配方施肥、養分分區宏觀調控和農業面源污染關鍵源區識別提供理論指導。

1 材料與方法

1.1 研究區概況

巢湖是我國五大淡水湖之一，水域面積780 km2，位于安徽省合肥市境內，屬北亞熱帶濕潤季風氣候。研究區域位于30°57'~32°16'N，116°40'~117°58'E，占地面積共8 607 km2（包含巢湖水域），耕地面積共3 556.7 km2（巢湖市700.2 km2、肥東縣1 054.2 km2、肥西縣674.9 km2、廬江縣1 105.4 km2、包河區22.0 km2），其中水田占比最大，為3 005.6 km2，占耕地面積的84.51%。研究區域為江淮丘陵地區，地形起伏變化明顯，土壤母質主要為下蜀系黃土，土壤類型以水稻土為主。

1.2 樣品采集與檢測

基于土地利用方式、行政區劃和耕地面積劃分采樣單元，并結合最新的GF 衛星影像數據和谷歌影像數據進行室內點位勘誤和糾正，最后通過點位實地踏勘確定采樣點位。研究區域共布點8 073 個（圖1），其中巢湖市1 632 個、肥東縣2 258 個、肥西縣1 555個、廬江縣2 295 個、包河區333 個。研究區域平均每0.465 km2耕地1個點位。

圖1 環巢湖地區采樣點分布圖Figure 1 Distribution of sampling points around Chaohu Lake

采樣時間為2021 年9 月至11 月，采樣方法按照《測土配方施肥技術規程》（NY/T 2911—2016），取表層（0~30 cm）樣5 點混合，多余樣品用四分法舍棄，最終取樣量在1~2 kg。對所有土壤樣品檢測pH、有機質、全氮、堿解氮、有效磷、速效鉀和緩效鉀，具體指標檢測規程和量綱如表1所示。

表1 土壤養分指標檢測規程Table 1 Methods for determination of soil nutrients

1.3 研究方法

1.3.1 數據預處理

為確保數據分析的準確性，需要對檢測數據進行剔除異常值和正態分布檢驗，由于該研究樣本量較大，剔除異常值采用3 倍標準差法（平均值±3 倍標準差）[11]，異常值用3 倍標準差范圍內的最大最小值替代，從而保證樣點數不變，并采用頻率直方圖的偏度和峰度來檢驗數據是否符合正態分布。

1.3.2 地統計分析

地統計方法是基于半變異函數來研究土壤元素空間變異特征的常用方法，能夠揭示隨機性和結構性因素對區域化變量的影響程度[12]。半變異函數公式如下：

式中：γ（h）為空間距離為h 時的半方差；N（h）為空間距離為h 的所有樣點的點對數；Z（xi+h）和Z（xi）分別為空間位置在xi+h和xi時的實測值。

半變異函數的參數中，塊金值（C0）代表人為活動和采樣檢測誤差等因素引起的隨機變異，基臺值（C0+C）為總空間變異程度，塊基比[C0/（C0+C）]為隨機變異在總空間變異中的占比，塊基比的大小能夠衡量空間自相關程度，其變程為空間自相關的理論距離。

1.3.3 空間自相關分析

空間自相關為空間上與相鄰區域之間觀測值的相互依賴性[13]，全局Moran's I 指數（IN）能夠反映整個研究區域的空間自相關性大小，其計算公式如下：

式中：N 為樣本數；xi和xj分別代表空間上在i處和j處的實測值；Wij代表空間權重；xˉ為樣本平均值。

全局Moran's I 指數在-1 到1 之間，大于0 表示正相關，小于0 表示負相關，其絕對值越大表明其自相關性越明顯，全局Moran's I 指數等于0 則表示不相關，空間分布為隨機分布。通過Z得分可以檢驗其顯著性，在0.01置信水平下，當|Z|≥2.58時表示為極顯著空間自相關。

1.3.4 空間聚類分析

空間聚類分析是基于局部Moran's I指數（Ii）來反映局部空間的自相關性，其計算公式如下：

局部Moran's I 最終將研究區域劃分為5 種情況，即①不顯著型：無顯著空間相關性，呈隨機分布；②高-高型：高值聚集區域；③高-低型：低值圍繞高值區域；④低-高型：高值圍繞低值區域；⑤低-低型：低值聚集區域。通過LISA 聚類圖直接反映出局部區域土壤屬性的關聯性[14]。

2 結果與分析

2.1 環巢湖地區土壤養分含量統計特征

對采樣點土壤檢測結果剔除異常值后的描述性統計結果及正態分布檢驗結果如表2 所示。土壤pH、有機質、全氮、堿解氮、有效磷、速效鉀和緩效鉀的均值分別為6.07、19.97 g·kg-1、1.16 g·kg-1、112.95 mg·kg-1、18.21 mg·kg-1、130.88 mg·kg-1和378.58 mg·kg-1，按照安徽省耕地質量監測指標分級標準，研究區土壤pH 為弱酸性，有機質、全氮、有效磷和速效鉀為中等水平，緩效鉀較低。根據變異系數小于10%為弱變異、10%~90%為中等變異、大于90%為強變異的劃分標準[15]，有效磷為強度變異，其變異系數達到了175.16%，其他土壤養分指標都為中等變異，其中pH變異系數最小，為13.60%。正態分布檢驗結果顯示pH、有機質、全氮和堿解氮符合正態分布，有效磷、速效鉀和緩效鉀經對數轉換后符合正態分布。

表2 環巢湖地區土壤養分描述性統計特征Table 2 Descriptive statistical characteristics of soil nutrients around Chaohu Lake

2.2 環巢湖地區土壤養分空間分布特征

2.2.1 地統計擬合模型的選擇

計算并選擇不同擬合模型下的半變異函數參數，是進行地統計分析和克里金插值的前提，表3 是不同擬合模型下的交叉驗證預測誤差。按照預測誤差的平均值和標準平均值接近0、均方根最小、標準均方根接近1、平均標準誤差接近均方根[11]，選擇各土壤養分的最佳擬合模型。結果表明土壤pH、全氮、有效磷、速效鉀和緩效鉀的最佳擬合模型為指數模型，有機質和堿解氮的最佳擬合模型為高斯模型。

2.2.2 土壤養分空間分布格局

在最佳擬合模型的半變異函數模型及參數的基礎上，使用ArcMap10.6 軟件通過普通克里金插值法制作環巢湖地區土壤各養分分布圖（圖2）。綜合7項土壤養分來看，環巢湖地區土壤養分都存在明顯的南北差異，具有不同程度的破碎化分布，空間變異特征明顯。表現在土壤pH 西南部偏酸性，東北部總體偏堿性，但在巢湖水面至東北方向有一條pH 低值條帶狀分布。土壤有機質、全氮和堿解氮的總體分布較為一致，由東南向西北方向其養分含量逐漸下降，并且在最北部有一塊明顯低值的三角形區域。有效磷小范圍分布破碎化較為明顯，整體上西部和北部地區偏低，但在西部和北部中間有一片顯著高值區域。速效鉀和緩效鉀都呈現為中北部顯著高于其他區域，從巢湖水面向東北方向有一條明顯的高低值分界線。

圖2 環巢湖地區土壤養分空間分布Figure 2 Spatial distribution map of soil nutrients around Chaohu Lake

2.3 環巢湖地區土壤養分空間變異特征

2.3.1 土壤養分地統計分析

在地統計分析中，步長大小的選擇十分重要，步長過大可能會掩蓋短程的自相關性，步長過小會導致落在計算范圍內的點對太少，從而影響研究結果[16]。本研究通過計算采樣點與最相鄰點之間的平均距離（451.04 m），采用略低于該值的步長（450 m）和系統默認的12 個步長數，從而確保所有的步長都會有若干個點對。環巢湖地區土壤養分擬合模型半變異函數參數如表4 所示，其中塊基比代表了隨機因素引起的空間變異在總變異中的占比，用來表示空間自相關程度，塊基比小于25%表示強烈的空間自相關性，25%~75%為中等的空間自相關性，大于75%則空間自相關性較弱[17]。環巢湖地區土壤pH 和緩效鉀為中等空間自相關性，表明空間變異由隨機性因素（施肥、種植習慣等人為活動）和結構性因素（母質、地形等自然條件）共同影響。其他養分指標（有機質、全氮、堿解氮、有效磷、速效鉀）空間自相關性均較弱，表明其變異主要受隨機性因素影響。變程代表區域變量的空間變異范圍，7 項養分指標的變程依次為堿解氮>有效磷>全氮>有機質>緩效鉀>速效鉀>pH，堿解氮的變程為4 978 m，說明在4 978 m 的范圍內其空間自相關性較弱。pH 的變程最小（1 968 m），說明環境因素在小范圍上影響著pH 的空間異質性，較其他土壤養分空間分布更加破碎化[18]。

表4 環巢湖地區土壤養分半變異函數參數與全局Moran′s I指數Table 4 Semivariogram parameters and global Moran′s I index for soil nutrients of around Chaohu Lake

2.3.2 土壤養分空間自相關分析

全局Moran's I 指數反映了區域變量在整體上的聚集狀態，其絕對值越大，說明空間自相關性越明顯。環巢湖地區土壤各養分均為極顯著空間正相關（Mo?ran's I>0，P<0.01，Z>2.58），其全局Moran's I指數依次為有效磷>速效鉀>pH>緩效鉀>全氮>有機質>堿解氮，有效磷的全局Moran's I值為0.711，說明采樣點有效磷含量與鄰近采樣點接近，其空間聚類更為明顯。有機質、全氮和堿解氮的全局Moran's I 值較小，表明其在空間上的分布更為分散。

2.3.3 土壤養分空間聚類分析

基于普通克里金插值結果，以行政村為評價單元，將行政村內插值結果的平均值空間連接到評價單元，并進行局部Moran's I指數計算，形成LISA 聚類圖（圖3）。pH高值聚集區域主要分布在北部和東部，且伴有少部分的低-高型異常區域，低值聚集區域主要分布在南部。有機質、全氮和堿解氮的聚類特征較為分散，其低值聚集區域都主要分布在中北部，高值聚集區域較為破碎化，在南部和中部均有不同程度的分布，同時高-低型和低-高型異常區域在高、低值聚集區域周邊都有零散分布。有效磷的空間聚集特征最為明顯，低值聚集區域在北部和西部，北部與西部中間有片明顯的高值聚集區域，且異常值區域較少。速效鉀和緩效鉀的高值聚集區域分布在中北部，低值聚集區域分布在中部和南部，且在高低值聚集區周邊伴隨少量的異常值區域。

圖3 環巢湖地區土壤養分LISA聚類分布圖Figure 3 LISA map of soil nutrients around Chaohu Lake

3 討論

3.1 養分變異特征分析與分區管理建議

3.1.1 土壤pH

環巢湖地區土壤pH 的塊基比（73.92%）相比較其他養分較小，其pH 的變程（1 968 m）和變異系數（13.60%）在各養分中最低，這與董佳琦等[19]的研究結果相似，表明pH的空間連續性較小，在小范圍內具有中等程度的空間自相關性，在土壤pH 分布圖上表現為小塊斑狀聚集分布。研究區整體上自東北至西南方向土壤逐漸偏酸性，這可能與其有機質南高北低的大致分布有關[20]，腐殖質分解的腐植酸會降低土壤pH[21]，針對偏酸性地區應選擇適宜的種植模式、保持水土、減少單質化肥的投入[22-23]，并且施用石灰能夠調節土壤pH 從而減少重金屬的危害[24]。地形地貌能夠通過影響成土母質、土壤類型、土地利用方式及水熱等條件，從而影響土壤屬性，研究區自巢湖水面至東北方向有一條帶狀分布的土壤pH 低值聚集區，此處為大別山余脈，其山地分布與土壤pH 低值聚集區呈現高度重合，這主要是由于山地地區土壤淋溶作用較強，鹽基離子易被氫離子取代成為鹽基不飽和土壤，從而導致土壤偏酸性[25]，對于山地酸性土壤地區可適當選擇種植果樹與煙葉等，并增施有機肥來提高土壤肥力。

3.1.2 土壤有機質、全氮和堿解氮

土壤有機質、全氮和堿解氮的相關性較高，空間分布特征較為相似，其塊基比（93.99%、91.39%、90.92%）和變程（3 874、4 467、4 978 m）較高，而全局Moran's I 指數（0.153、0.155、0.107）較小，表明在較大范圍內人為活動等隨機性因素主導其空間變異，且空間分布較為分散，在其養分分布圖中呈現為小塊斑狀零散分布，這是由于土壤有機質受基施有機肥的影響較大[26]，小農戶經營的現狀致使種植模式難以統一，從而導致了在一定范圍內土壤有機質、全氮和堿解氮分布破碎化。土壤有機質、全氮和堿解氮的變異系數（37.89%、36.91%、44.97%）相對較小，說明整體離散程度較低，具有破碎化嚴重又向均質化發展的分布特征。因此，在測土配方工作中應適當提高樣點密度，從而更精確地掌握該養分分布特征，促進化肥減量增效。在土壤養分指標中，有機質和全氮含量對作物產量影響較大，其與土地利用方式和有機肥的施用相關性較高[27]，針對有機質和全氮偏低的地區實行有機無機配施、作物輪作[28]、秸稈還田和間套種植模式[27]，能夠有效改善土壤理化性質并提高有機質與有機碳等含量。

3.1.3 土壤有效磷

有效磷的塊基比（85.68%）和變程（4 878 m）與有機質（和全氮、堿解氮）較為接近，說明有效磷在大范圍內的強度變異主要由人為活動等隨機性因素導致[26]。其變異系數（175.16%）和全局Moran's I 指數（0.711）顯著高于其他養分，且空間自相關性極為顯著（Z>2.58），表現為高低值聚類特征明顯，整體離散程度也最大，綜合說明了有效磷具有高度積累的現象。結合LISA 圖了解到有效磷在西部與北部有明顯的大片低值聚集區域，但在兩地中間卻有一片高值聚集區域，這可能是中間區域為合肥市城區，人為活動較為集中，導致有效磷富集嚴重。土壤中的磷主要以吸附態和磷酸鹽為主，其移動性較弱，極易被固定，當季利用率較低[29]。針對有效磷集中偏高的市區，應當尋找導致其含量普遍偏高的具體原因，采取“源頭控制-過程攔截-末端治理”的防控思路并結合政策調控，從而避免造成農業面源污染[30]；而對有效磷偏低的地區，搭配有機肥與磷肥配施能夠提高磷的有效性，并且施用磷肥能夠明顯提高有效磷含量[31]，同時有機質也在一定程度上能夠影響有效磷的含量[32]。

3.1.4 土壤速效鉀和緩效鉀

速效鉀和緩效鉀在空間分布上保持著高度統一，但速效鉀的塊基比（81.77%）明顯高于緩效鉀（63.20%），說明速效鉀的空間變異由人為活動等隨機性因素的影響相對較大，而緩效鉀更多的由成土母質等結構性因素影響，這是由于土壤中的速效鉀為交換性鉀和水溶性鉀，而緩效鉀為黏土礦物固定的鉀，因此與成土母質的關系更為密切[33]。同時，速效鉀的全局Moran's I 指數（0.397）也高于緩效鉀（0.213），體現為速效鉀的空間自相關性更為明顯，更容易呈現聚集性分布，這是由于當作物吸收而導致速效鉀不足時，緩效鉀會迅速轉化為速效鉀，從而維持速效鉀的含量。與土壤pH 分布特征類似，速效鉀和緩效鉀的含量在南部丘陵地區與東北方向的條帶狀山地區偏低，這是由于海拔較高的低山丘陵區受水流侵蝕作用較強，較土壤氮、磷元素，鉀元素更容易流失[13]。作物秸稈中含有大量鉀素，對缺鉀地區促進秸稈規范化還田，注重施用鉀肥能夠有效提高土壤鉀含量。

3.2 研究創新點與改進思路

本研究結果中，土壤養分空間自相關性總體上較弱，其空間變異更多的是由人為活動等隨機性因素導致，且主導了土壤養分變異朝均質化方向發展[19]。在地統計分析結果中，土壤pH更容易受到母質、地形等結構性因素的影響，這與其他地統計學的研究結論一致[8，34]，并且土壤pH、速效鉀和緩效鉀的分布與地形海拔相關性較高，有許多學者已經證實了這一點[13，35]；有效磷的變異系數較其他養分更大，這與牛文鵬等[36]在多養分中的對比結果一致，有效磷的全局Moran's I 指數最大也證實了其更容易富集的特點。有機質、全氮和堿解氮的塊基比最大，這與姜霓雯等[13]的研究中有機質與全氮塊基比最大的結果類似，說明了其在土壤養分中更容易受到人為活動等隨機性因素影響。從研究數據上看，本研究數據與其他相關研究偏差較大，但養分之間差異的總體規律較為相似，這可能與本研究區布設樣點數較多有關。這是由于區域變量的空間變異特征結果隨著采樣點尺度變化而變化[16]，當采樣點的尺度明顯大于母質、地形等自然因素分布的尺度時，地統計分析結果更偏向于由結構性因素主導其空間變異[36]，當采樣點尺度較小時結果則相反，在部分高密度采樣的研究結果中已經證實了這一點[37]。

本研究基于環巢湖地區大范圍尺度高密度隨機布點采樣的檢測數據，使用半變異函數參數的塊基比和變程分別說明其養分空間自相關程度和作用范圍，用全局Moran's I 指數反映其空間自相關性的大小及顯著性，并將研究結果與養分統計特征的變異系數呼應，最后通過克里金插值和局部Moran's I 指數制作養分分布圖和LISA 聚類圖，系統地呈現了研究區多種土壤養分的空間變異特征及內在聯系。雖然利用地統計學和Moran's I 指數研究土壤養分空間變異特征的研究較多，但同時使用兩種方法交叉進行結果論證的研究卻較少，同時本研究涉及到的土壤養分較為全面，進一步對比分析更能為其他土壤養分空間變異相關研究提供可靠的數據參考。本研究詳細介紹了環巢湖地區土壤養分的空間變異特征，但未深入探討養分與種植模式、土壤質地、坡度坡向、海拔和植被類型等多種驅動力之間的關系，隨著遙感技術的發展，下一步研究可以通過遙感影像提取或預測這些影響因素，從而更好地了解土壤養分空間分布的關鍵因素以及發展趨勢。

4 結論

（1）環巢湖地區土壤pH、有機質、全氮、堿解氮、有效磷、速效鉀和緩效鉀的均值分別為6.07、19.97 g·kg-1、1.16 g·kg-1、112.95 mg·kg-1、18.21 mg·kg-1、130.88 mg·kg-1和378.58 mg·kg-1，除有效磷為強度變異外，其他養分均為中等程度變異。

（2）環巢湖地區土壤pH、全氮、有效磷、速效鉀和緩效鉀的半變異函數最佳擬合模型為指數模型，有機質和堿解氮的最佳擬合模型為高斯模型。pH 和速效鉀為中等空間自相關性，其他養分空間自相關性均較弱，表明土壤養分空間變異主要受人為活動等隨機性因素影響。

（3）環巢湖地區土壤各養分均為極顯著空間正相關，其空間相關性大小依次為有效磷>速效鉀>pH>緩效鉀>全氮>有機質>堿解氮，且有效磷的空間自相關性顯著大于其他養分。

（4）環巢湖地區土壤除有效磷外的其他養分空間分布較為分散，且不同養分空間分布規律不同，總體上看，南北方向各養分含量逐漸上升或下降，其養分高、低值聚集區域也有明顯的南北分布。而有效磷聚集特征明顯，表現為西部至北部兩端低值聚集，中間區域高值聚集。