基于GIS的濱海鹽漬化農田土壤空間變異及其分區管理

朱昌達, 高明秀,*, 王文倩, 李俊翰, 姚 宇, 周娜娜

1 山東農業大學資源與環境學院, 泰安 271018 2 土肥資源高效利用國家工程實驗室, 泰安 271018

精準農業是現代農業發展的方向,長期以來備受關注。全面分析掌握農田土壤特征,精準劃分管理區,是采取差異化管理措施,實現精準農業的基礎[1-2]。近年來,一些學者針對精準管理分區方法展開研究,探討了運用決策樹、系統聚類、K均值聚類(灰色聚類)和模糊聚類等方法[3-6]劃分管理分區的可行性,其中K均值聚類和模糊c-均值聚類(FCM,Fuzzy c-means algorithm)[7-8]應用較多。由于很多分類對象間沒有明確的分類界線、存在亦此亦彼的表現,因此相比而言,模糊c-均值聚類作為一種新的空間連續性聚類算法,克服了K均值聚類只能將每個對象嚴格地分類到某一聚類中心、與其他聚類中心的空間相關關系得不到體現的缺陷[9],對于復雜的土壤管理問題更為適宜。基于FCM劃分管理分區的關鍵在于聚類分區數目c和模糊加權指數φ的確定,已有研究主要采取經驗取值法、c-φ多次組合最優法、指數檢驗法等[10-12],對分區結果進行檢驗主要分為基于內部變量和外部變量兩種方法[7,13]。

綜觀現有農田精準管理分區研究,小尺度的試驗田研究較多、大中尺度的農田研究少,土壤條件較好的內陸田塊試驗較多、濱海鹽漬化區域缺乏涉及；只針對單一土壤屬性的研究較多,綜合考慮土壤多個屬性的研究較少；從分區結果檢驗來看,基于單方面變量的檢驗都存在較明顯的誤差[14]。因此,本文擬以山東省無棣縣農田為例,針對濱海鹽漬化農田鹽堿瘠薄、土壤屬性空間變異大、粗放管理效益低的現實,在縣域尺度上,研究基于多個土壤屬性變量的農田管理分區精準劃分方法和分區精度校驗方法,以提高分區結果的合理性和科學性,為濱海鹽漬化地區采取差異化管理措施,提升鹽漬化土地利用水平提供科學依據[15]。

1 材料與方法

1.1 研究區概況

無棣縣(37°41′—38°18′N,117°31′—118°12′E)地處山東省最北部,總面積1963 km2,東北部瀕臨渤海,屬北溫帶東亞季風區域大陸性氣候,年均降水570.1 mm,80%集中在夏季,秋冬春三季干旱,年均蒸發1285.5 mm,蒸降比2.3。近年來,該縣作為濱海低平原區新增耕地潛力和增產潛力較大的縣區之一,置身新增千億斤糧食生產能力規劃、鹽堿地治理、“渤海糧倉”科技示范工程等國家戰略疊加區,成為農業開發熱點區域[16]。但是,無棣縣淡水匱乏,地下水位淺、礦化度高,土壤瘠薄鹽堿問題突出[17],農業生產效率和效益較低,如何通過精準分區、差異化管理提高農田利用水平和生產效益成為亟需解決的問題。由于縣境東北部主要為鹽田、海水養殖區和荒草地,因此本文研究區僅涉及農田分布區,總面積約為756 km2。

1.2 數據來源

在土地利用現狀圖上以2 km×2 km網格覆蓋研究區農田(避開建設用地、鹽田、灘涂),面積不足者并入相臨網格,每個網格中心設計1個調查樣點(圖1),野外調查時根據現場情況確定實際采樣地點,使用手持式GPS記錄坐標,并用五點取樣法采取對應的土壤樣品1 kg左右裝袋密封,帶回實驗室化驗。野外調查在2016年5月24—28日進行,共采土樣193個。土樣經自然風干、研磨、過篩,化驗分析[18]獲得193組土壤有機質、有效氮、有效磷、速效鉀、含鹽量和pH等屬性數據。

圖1 無棣縣網格劃分及采樣點位分布圖 Fig.1 Grid division and sampling point distribution map in Wudi County

1.3 研究方法

1.3.1土壤屬性特征與變異性分析

運用SPSS 20.0進行土壤屬性描述性統計分析,并通過Person-雙側檢驗得到相關系數矩陣,分析土壤各屬性變量的變異性及相關性。在ArcGIS 10.2中運用地統計學方法,分析土壤屬性的空間分布特征。土壤養分等級參照全國第二次土壤普查分級標準劃分為極低、低、中、豐、高5級,含量區間從低至高依次為,有機質(g/kg)：(0,6)、[6,10)、[10,20)、[20,30)、[30,∞)；有效氮(mg/kg)：(0,30)、[30,60)、[60,90)、[90,120)、[120,∞)；有效磷(mg/kg)：(0,5)、[5,15)、[15,30)、[30,40)、[40,∞)；速效鉀(mg/kg)：(0,40)、[40,75)、[75,120)、[120,150)、[150,∞)；含鹽量(g/kg)結合作物耐鹽能力分為脫鹽化、輕度、中度、重度和極重度5級：(0,1)、[1,2)、[2,3)、[3,4)、[4,∞)；pH分為酸性、適宜、中等堿性、強堿性和極強堿性5級：(0,6.5)、[6.5,7.5)、[7.5,8.5)、[8.5,9)、[9,∞)[17]。

1.3.2模糊均值聚類

模糊c-均值聚類(FCM)是一種非監督的連續性動態聚類算法,定義一個n×p的土壤屬性多源數據集(n=193,p= 6),最常用的目標函數為：

(1)

聚類有效性檢驗通過模糊性能指數(FPI,Fuzziness performance index)和歸一化分類熵(NCE,Normalized classification entropy)評價。FPI和NCE的取值范圍在0—1之間,取值越接近0,分類界線越明顯[19],反之分類界線越模糊。公式為：

(2)

(3)

式中,H(U;c)是分類熵的函數表達式,a可以取任意正整數。

1.3.3管理分區劃分及可視化

運用MAPGIS 6.7進行研究區底圖及采樣點繪制,并將得到的隸屬度賦屬性到各采樣點中,實現管理分區的可視化表達。由于在MAPGIS中只能依據最大隸屬度原則進行分類,不能體現模糊連續分類在空間預測中的重要作用,因此運用ArcGIS地統計分析模塊對土壤屬性模糊隸屬度進行插值分析,實現其空間連續性表達。

1.3.4分區精度檢驗

為了評價基于模糊聚類分析進行管理分區的精度,確定其變異程度是否滿足管理分區的要求,對不同管理分區內的土壤屬性進行變異性分析,并用最小極差法(LSR,Least significant range method)對各分區結果進行差異顯著性檢驗,校驗管理分區結果的精度。

2 結果與分析

2.1 土壤屬性特征及其變異性

在SPSS 20.0支持下,采用鄰域法對獲取的193組土壤屬性數據進行特異值處理[20],并進行描述性統計分析(表1)。由表1可見,研究區農田土壤有機質含量均值19.253 g/kg,處于中等水平；有效氮均值59.349 mg/kg,處于較低水平；有效磷均值29.527 mg/kg,處于中等水平；速效鉀均值179.676 mg/kg,含量較高；含鹽量均值1.204 g/kg,總體為輕度,這與采樣時已進入夏季、土壤水分含量較高有關；pH均值8.184,為中度堿性,土壤鹽漬化程度總體呈輕至中度水平。有機質、有效氮、有效磷、速效鉀和含鹽量均呈中等變異性(CV< 10%為弱變異性,10% ≤ CV< 100%為中等變異性,100% ≤ CV為強變異性),pH值總體較高但變異性較弱,是否作為管理分區的依據有待于分析與其他土壤屬性的相關性。綜合來看,研究區土壤屬性總體呈中等變異,有必要進行分區管理。

表1 土壤屬性描述性統計Table 1 Descriptive statistics of soil properties

2.2 土壤屬性相關性分析

在SPSS 20.0中對土壤屬性變量進行相關性分析,得到相關系數矩陣(表2)。由表2可得,有機質、有效氮、有效磷、速效鉀和pH值等5個指標的相關性較高。含鹽量雖相關性較低,但由于研究區位于環渤海低平原區,地下水埋深淺且礦化度高,土壤鹽漬化風險大,有必要將含鹽量作為分區評價指標之一。因此,確定將上述6個土壤屬性作為管理分區的評價指標。

表2 土壤各屬性的相關系數矩陣Table 2 Correlation coefficient matrix of soil properties

2.3 土壤屬性空間變異特征

在SPSS 20.0中用單樣本K-S(Kolmogorov-Smironov)工具進行正態分布檢驗,選擇默認顯著性為0.05、置信區間為95%,若數據的漸進顯著性P值>0.05,則假設成立,符合正態分布(表3)。由表3可見,有機質、有效氮、pH的K-S值均大于0.05,為正態分布；有效磷、速效鉀和含鹽量為偏態分布,經對數變換后K-S值大于0.05,達到顯著水平,呈正態分布,可以進行克里格插值分析。

表3 土壤屬性正態分布檢驗Table 3 Normal distribution test of soil properties

運用ArcGIS 10.2地統計工具進行克里格插值,采用C-V交叉檢驗法,通過多次試驗選擇最優擬合模型[21]。最終確定有機質、速效鉀的擬合模型均為指數模型,有效磷、有效氮、含鹽量和pH的擬合模型為球面模型。此時,土壤屬性各擬合模型的標準平均值最接近于0,標準均方根最接近于1(表4)。

表4 土壤屬性變異函數模型及相關參數Table 4 Variation function model and related parameters of soil properties

一般來說,塊金效應(Nugget / Sill)的取值小于25% 時,具有強空間相關性；在25%—75%之間時,具有中等空間相關性；大于75%時,具有弱空間相關性。由表4可見,土壤有機質、有效氮和有效磷的塊金效應值分別為51.9%、67.4%、67.0%,說明這三種土壤屬性的空間相關性為中等水平,空間變異特征受隨機因素和結構因素的共同影響。速效鉀、含鹽量的塊金效應值分別為20.8%、22.4%,都小于25%,說明含鹽量和速效鉀具有強空間相關性,其變異特征受隨機因素影響較小,主要由結構性因素引起。pH的塊金效應值為0,說明pH具有完全的空間自相關性,其空間變異特征完全由結構性因素引起。

運用克里格法確定土壤屬性最佳擬合模型后,生成土壤屬性的空間插值圖(圖2),以對土壤屬性空間變異規律進行最優預測,實現土壤屬性空間變異的數量化及可視化。由圖2可見,有機質高值區主要分布在研究區西南和西北部,呈塊狀分布。有效氮和有效磷高值區的分布趨勢較為一致,中部和南部分散分布,含量最高或最低處呈小塊狀或點狀分布。速效鉀的分布具有明顯的地帶性特征,高值區主要分布在中東部地區,呈東高西低現象。含鹽量在中部地區較低,且分布范圍廣,含量約為1.0—2.0 g/kg,屬于輕度鹽漬化,少數點狀分布地區含鹽量在2.0—4.0 g/kg之間,屬于中重度鹽漬化。pH值呈塊狀分布,值高(8.5—9.0)的地方主要分布在中東部地區。

圖2 不同土壤屬性的空間預測插值圖Fig.2 Spatial prediction interpolation of different soil properties

2.4 管理分區劃分及其可視化

在MATLAB R2016a中通過“z-score”數據標準化處理,得到的n×p(n=193,p=6)的土壤屬性模糊數據集。調用FCM函數,取冪指數為1.6,最大迭代次數為200,收斂閾值0.000001,進行模糊聚類分析。為了確定最佳分類數c,選取分類數2—8分別進行模糊聚類分析,并根據模糊性能指數(FPI)和歸一化分類熵(NCE)進行評價(圖3),兩者的值越小,聚類效果越明顯。由圖3可以看出,當分類數c取值為3時,FPI和NCE的值最小,說明此時同一分類中的數據差異性最小,不同分區間的差異性最大,確定本研究區的最佳分類數為3個。

圖3 不同分類數的FPI和NCE值 Fig.3 FPI and NCE values of different classification numbers FPI: 模糊性能指數 Fuzziness performance index; NCE: 歸一化分類熵 Normalized classification entropy

由模糊聚類分析得出的隸屬度矩陣,依據最大隸屬度原則將采樣點分為3類,在MAPGIS 6.7中根據各采樣點的最大隸屬度劃分管理分區,實現管理分區的空間可視化(圖4)。由于最大隸屬度原則是以網格為單位進行分類,不能完全體現模糊分類的空間預測性。本文根據分類后各采樣點的模糊隸屬度矩陣,在ArcGIS 10.2中利用地統計模塊進行插值分析,實現土壤屬性模糊隸屬度的空間連續性表達(圖5)。

圖4 基于模糊聚類分析的管理分區圖 Fig.4 Management zones map based on fuzzy clustering analysis

圖5 各分區土壤模糊隸屬度值空間預測分布圖Fig.5 Spatial prediction zones map of soil fuzzy membership value in different zones

對比圖5和圖2可得,分區Ⅲ主要分布于有機質含量豐等級、有效氮含量中和豐等級、有效磷含量高等級區域；分區Ⅰ的中西部主要分布于有機質含量中等級區域,東南部地區主要分布于有效氮含量豐等級區域；分區Ⅱ主要分布于速效鉀含量高等級區域,其他養分含量較低,可以考慮多施鉀肥除外的其他肥料。且鹽漬化嚴重的北部地區和pH強堿性區域主要為分區Ⅰ和分區Ⅱ,與分區Ⅲ土壤的擬合度較低。因此,管理分區和土壤各屬性的空間分布具有很強的空間相關性和空間擬合度,各分區內土壤屬性向均一化方向發展。對比圖4和圖5可得,各分區土壤模糊隸屬度空間預測圖與最大隸屬度原則劃分的管理分區圖在空間分布上基本一致,各分區間的隸屬度較為明顯,交叉重疊程度低,說明以最大隸屬度原則劃分的管理分區具有空間預測的有效性。

統計各分區土壤屬性平均值(表5)可見,分區Ⅲ中的有機質、有效氮、有效磷含量均最高,速效鉀含量較高,土壤養分水平總體最高,且含鹽量、pH值較低；分區Ⅱ中的速效鉀含量最高,有機質、有效氮的含量在分區Ⅲ和分區Ⅰ之間,有效磷比分區Ⅰ略低,土壤養分水平總體居中；分區Ⅰ的土壤養分水平總體最低,含鹽量和pH值也最高。

2.5 管理分區精度校驗

為了對基于模糊聚類分析進行管理分區的結果進行精度驗證,對不同分區內土壤屬性進行統計分析,得到基于變異系數和最小極差法(LSR)的各分區差異顯著性檢驗結果(表5)。

對比表1和表5可見,各分區的土壤屬性變異系數較全研究區均有所下降。其中,有機質變異系數由25.0%下降為14.7%—23.9%,有效氮變異系數由30.0%下降為17.9%—30.0%,有效磷變異系數由52.3%下降為33.0%—51.5%,速效鉀變異系數由32.8%下降為18.2%—29.3%。pH值變異系數由3.0%下降為2.7%—3.0%,變化較小,均屬于弱變異性,分區結果對pH影響較小；含鹽量變異系數由48.4%變為41.0%—54.1%,分區Ⅲ變異系數下降,分區Ⅰ和分區Ⅱ均略有升高,說明含鹽量分布不均勻。

表5 各分區土壤屬性描述性統計和LSR檢驗結果Table 5 Descriptive statistics and LSR test results of soil properties in different zones

LSR檢驗結果表明(表5),土壤有機質、有效氮、速效鉀在3個分區均達到極顯著差異(P<0.01)或顯著差異(P<0.05),有效磷含量在分區Ⅲ和分區Ⅱ之間具有顯著差異(P<0.01),含鹽量在分區Ⅱ和分區Ⅰ、分區Ⅲ間具有極顯著差異(P<0.05),pH值在分區Ⅲ和分區Ⅰ、分區Ⅱ之間具有極顯著差異(P<0.01)。

總體而言,分區內的土壤屬性變異程度較分區前全區的變異程度有所下降,分區間的差異顯著性明顯,分區結果具有可行性。其中分區Ⅲ的變異程度下降最明顯,各采樣點土壤屬性隸屬度最明確,分區Ⅰ下降程度最低,分區Ⅱ居中。基于變異系數和LSR的檢驗結果均表明,管理分區具有較高的精度。各分區內部的土壤養分變異系數由25%—52.3%減小為14.7%—51.5%,各分區內部的空間變異性較全研究區減小,且各分區間具有顯著性差異,分區結果可以作為農田分區調控的統一作業單元。根據圖5統計各分區面積,Ⅰ、Ⅱ、Ⅲ管理區面積分別為2.56萬hm2、1.76萬hm2、3.24萬hm2。

3 結論

該文以無棣縣農田為研究區,采用網格法結合土地利用現狀定點野外采樣、室內化驗分析獲取土壤屬性數據,運用ArcGIS 10.2地統計方法分析土壤屬性的空間變異特征；在MATLAB R2016a中采用模糊c-均值聚類法(FCM)計算各樣點的模糊隸屬度,通過插值預測模糊隸屬度的空間分布,基于最大隸屬度原則進行分區,并根據模糊隸屬度分布圖分析分區結果的有效性；通過變異性分析和最小極差法(LSR)差異顯著性檢驗,對分區結果進行精度驗證。

研究結果表明：無棣縣農田土壤總體呈輕至中度鹽漬化,有效氮含量偏低,有機質、有效磷含量中等,速效鉀含量較高；有機質、有效氮、有效磷、速效鉀和含鹽量呈中等變異性(變異系數25%—52.3%),空間變異性較大,應分區調控；空間變異分析顯示,速效鉀、含鹽量和pH的塊金效應值都小于25%,受人為因素影響相對較小,主要受土壤質地、地下水礦化度等結構因素影響,有機質、有效氮和有效磷的塊金效應值在50%—75%之間,受耕作方式、施肥等人為因素影響較大；將全縣農田劃分為Ⅰ、Ⅱ、Ⅲ共3類管理區,估算面積分別為2.56萬hm2、1.76萬hm2、3.24萬hm2；各分區土壤養分變異系數分別為23.9%—51.5%、15.9%—50.3%、14.7%—33.0%,檢驗結果表明各分區間差異顯著,而各分區內部變異性明顯低于全研究區。管理分區與土壤屬性空間分布特征具有較高的擬合度,表明分區結果可以作為差異化管理的作業單元。研究結果為各分區內部統一、不同分區間差異化管理提供了依據,研究有助于推進濱海鹽漬化農田精準化管理水平的提高。

由于分區后pH和含鹽量的變異系數變化不大,少數區域甚至出現含鹽量變異性增大的現象,考慮到當地受濱海鹽漬化影響嚴重,土壤水鹽狀況短時間難以改變,可在上述管理分區的基礎上,對變化不明顯或者變異性增大的區域進行細化管理,以更好地改善其水鹽狀況。