基于Moran′s I的菜地土壤屬性空間分布格局分析

王強，鄭夢蕾，葉治山，楊善蓮，馬友華

（1.安徽農(nóng)業(yè)大學(xué)資源與環(huán)境學(xué)院，合肥230036；2.安徽農(nóng)業(yè)大學(xué)新農(nóng)村發(fā)展研究院，合肥230036）

土壤養(yǎng)分累積遷移成為地表水和地下水的污染源，是全世界迫切解決的難題之一[1]。農(nóng)田面源污染中氮磷元素的遷移已導(dǎo)致我國大面積的湖泊和水庫富營養(yǎng)化[2]，受到了普遍關(guān)注，一些學(xué)者利用多元統(tǒng)計、地統(tǒng)計方法進行定性或定量研究[3]。一般認為自然來源和人為來源是農(nóng)田面源污染的兩大來源，前者指在成土過程中地質(zhì)高背景值帶來的富集作用，后者主要指農(nóng)業(yè)生產(chǎn)及農(nóng)村生活對土壤的作用[4-7]。蔬菜在居民生活中必不可少，我國已是世界上最大的蔬菜生產(chǎn)和消費國，蔬菜播種面積和產(chǎn)量均占世界的40%以上[8]。根據(jù)中國國家統(tǒng)計局數(shù)據(jù)及公開發(fā)表文獻，我國蔬菜面積已達到1 998.1 萬hm2，設(shè)施蔬菜面積370 萬hm2[9]，蔬菜露天栽培面積1 628.1 萬hm2。我國蔬菜集約化種植區(qū)大量施用肥料，導(dǎo)致土壤養(yǎng)分流失，引起土壤酸化、鹽漬化和地下水污染，使土壤生物活性降低、分解能力下降，造成了嚴重的經(jīng)濟及環(huán)境問題[10]。農(nóng)業(yè)景觀空間格局是農(nóng)業(yè)面源污染的主要影響因素之一[11]，菜地面源污染風(fēng)險空間分布格局識別與優(yōu)化可為科學(xué)防范和治理農(nóng)業(yè)面源污染提供決策依據(jù)[12-14]。農(nóng)業(yè)景觀空間格局主要研究景觀格局的時空變化特征[15]，通常采用遙感、地理信息科學(xué)和景觀格局分析方法，運用經(jīng)典回歸模型和空間滯后模型研究農(nóng)業(yè)景觀格局及其組分的變化，如傅伯杰[16]對黃土區(qū)農(nóng)業(yè)景觀空間格局的研究，Liu 等[17]運用多元線性回歸模型量化不同土地利用對流域中氮、磷濃度貢獻的研究。目前國內(nèi)外學(xué)者在農(nóng)業(yè)景觀研究領(lǐng)域已有相關(guān)研究，但主要局限于農(nóng)業(yè)景觀格局變化及其環(huán)境效應(yīng)、評價方法及其應(yīng)用等方面[18]?，F(xiàn)代農(nóng)業(yè)要求農(nóng)田集約化、田面平整化、田塊規(guī)則化，這促使農(nóng)田斑塊均質(zhì)，導(dǎo)致土壤生態(tài)系統(tǒng)自調(diào)節(jié)功能減弱，為探明生態(tài)風(fēng)險區(qū)來源，需要可追溯污染源的空間分析方法[19]。局部二元莫蘭指數(shù)分析可以幫助檢測空間異常值或熱點，是一種有效的區(qū)域變量探索性空間分析方法[20]。一些學(xué)者研究了我國菜地土壤養(yǎng)分的時空變化和環(huán)境質(zhì)量變化[10]。Valente 等[21]使用局部二元莫蘭指數(shù)進行空間相關(guān)性分析，結(jié)果表明土壤電導(dǎo)率與土壤未利用磷呈強相關(guān)性，可作為精準農(nóng)業(yè)管理的工具。Fu 等[22]研究發(fā)現(xiàn)，由于農(nóng)民的集約化管理，土壤磷在農(nóng)場周圍和交通路線呈現(xiàn)高-高（High-High）空間分布。然而以往研究中缺乏時空尺度下菜地種植對土壤屬性空間分布格局影響的研究，無法闡明研究區(qū)長時間尺度下菜地種植對土壤屬性值變化的影響機制。因此，本研究以安徽省肥東縣為實驗區(qū)，通過面向?qū)ο蠓诸惙椒ǐ@取2019 年農(nóng)業(yè)設(shè)施菜地空間分布位置，借助地面農(nóng)業(yè)調(diào)查得到研究區(qū)全部菜地空間分布位置，結(jié)合采集的375 個表層土壤樣品，運用空間插值方法模擬土壤屬性含量空間分布，運用莫蘭指數(shù)對菜地空間密度分布是否影響土壤屬性空間分布格局的問題進行研究，研究結(jié)果為農(nóng)業(yè)土壤質(zhì)量安全的管理與保護提供數(shù)據(jù)支撐，有助于基本農(nóng)田的長期管理和保護。

1 材料和方法

1.1 研究區(qū)概況

肥東縣位于安徽省中部，是巢湖流域面源污染優(yōu)先控制區(qū)[23]，面積22.12 萬hm2，耕地面積7.67 萬hm2，占總面積的35%，總?cè)丝?10萬，其中農(nóng)業(yè)人口97萬，是一個典型的農(nóng)業(yè)大縣。氣候?qū)俦眮啛釒Ъ撅L(fēng)氣候區(qū)，四季分明，雨量適中。菜地總面積6 159 hm2，其中露天菜地5 092.52 hm2，種植時間從1989 年至2013年，設(shè)施菜地1 066.48 hm2，種植時間從2003 年至2019年。

1.2 數(shù)據(jù)來源及技術(shù)路線

1.2.1 蔬菜種植類型獲取

遙感影像來源于2019年下載的Google earth17級數(shù)據(jù)，主要根據(jù)設(shè)施的形態(tài)特征如長度、寬度、空間分布等，基于面向?qū)ο蠓椒ㄌ崛∮跋裰械乃性O(shè)施菜地。運用谷歌地球中的不同歷史影像，對所有的設(shè)施菜地進行時空對比，將最早出現(xiàn)時間作為起始種植時間。設(shè)施菜地提取結(jié)果與2016 年全縣農(nóng)業(yè)調(diào)查總菜地空間分布進行空間相交，分別得到不同種植年限露天菜地與設(shè)施蔬菜空間分布，如圖1所示。

圖1 兩種菜地不同種植年限空間分布圖Figure 1 Spatial distribution of vegetable planting types

1.2.2 土壤采樣與分析

基于成土母質(zhì)、土壤類型劃分采樣單元并簡單隨機布點采樣，于2017 年7 月至12 月采集了375 個0～20 cm的表層土樣，包含24個露天菜地采樣點和10個設(shè)施菜地采樣點，如圖2 所示，其分析方法依據(jù)測土配方施肥技術(shù)規(guī)范（2011年修訂版）執(zhí)行。

圖2 采樣點分布圖Figure 2 Distribution of sampling points

1.3 地統(tǒng)計插值和核密度模擬

本文采用ARCMAP 10.2 軟件平臺進行土壤屬性的地統(tǒng)計學(xué)插值和菜地密度的模擬，地統(tǒng)計是依據(jù)協(xié)方差函數(shù)對隨機過程/隨機場進行空間建模和預(yù)測（插值）的一種回歸算法[24]。核密度估計是對離散柱狀圖連續(xù)替換，創(chuàng)建平滑曲線的估值方法[25]，把空間配置和輸入點數(shù)量計算的默認搜索半徑作為帶寬[26]，避免了稀疏數(shù)據(jù)集中的空間異常值。

1.4 空間自相關(guān)分析

莫蘭指數(shù)是計算空間自相關(guān)中最為經(jīng)典的方法[27]，一般分為全局型（Global spatial autocorrelation）和區(qū)域型（Local spatial autocorrelation）兩種[28]。全局型Moran′sI是基于統(tǒng)計學(xué)相關(guān)系數(shù)的共變量（Covariance）關(guān)系推算得來，共變量的大小程度即代表兩組數(shù)的相關(guān)性大小。全局型Moran′sI值的功能在于描述某現(xiàn)象的整體分布聚集狀況，5% 顯著水平下，Z（I）＞1.96時，表示研究范圍內(nèi)某現(xiàn)象的分布有顯著空間自相關(guān)性，Z（I）＜-1.96 時，表示研究范圍內(nèi)某現(xiàn)象的分布呈現(xiàn)負的空間自相關(guān)性。全局型的Moran′sI的公式如下：

式中：Wij是研究范圍內(nèi)每一個空間單元i與j（j={1，2，3，…，n}）區(qū)空間單元的空間相鄰權(quán)重矩陣，以1表示i與j相鄰，而以0 表示i與j不相鄰；xi為一組變量的每項值，xj為另一組變量的每項值，是其平均數(shù)。

局域型Moran′ sI是空間自相關(guān)局部檢驗統(tǒng)計量，用于識別空間聚類和空間異常值的位置。計算方法如下：

各變量定義與公式（1）類似。依照公式（2）計算出的Moran′sI值介于-1 到1 之間，大于0 為正相關(guān)，小于0 為負相關(guān)。值越大空間分布相關(guān)性越大，即空間聚集分布明顯，反之，值越小空間分布相關(guān)性越小，而當值趨于0 時，即代表此時空間分布呈現(xiàn)隨機分布的情形。采用GeoDa 軟件獲取全局和局部莫蘭指數(shù)空間自相關(guān)分類結(jié)果[27]。

2 結(jié)果與分析

2.1 菜地土壤屬性差異性分析

2.1.1 土壤屬性特征分析

根據(jù)前期調(diào)研，研究區(qū)青菜有機肥用量為11.25 t·hm-2，小白菜有機肥用量為7.5 t·hm-2，用量最多的肥料品牌為河南蓮鑫寶肥業(yè)有限公司的復(fù)合肥料（硫酸鉀型，總養(yǎng)分≥14%，有機質(zhì)≥45%）。土壤采樣數(shù)據(jù)的分析結(jié)果如表1 所示，土壤有效磷的變異系數(shù)最高，速效鉀、全氮、有機質(zhì)及pH 緊隨其后，分別達到29.34%、17.73%、14.07%、12.35%、4.67%。

如表2 所示，距離城鎮(zhèn)越近，有機質(zhì)、全氮含量和pH 相對越高，隨著距離的增加整體呈下降趨勢；有效磷含量和速效鉀含量則隨著距離的增加呈現(xiàn)先增加后下降的趨勢，在2.0～2.5 km處出現(xiàn)最高值。

如表3 所示，從均值可知隨著種植年限的增加，露天菜地土壤有機質(zhì)、全氮、有效磷、速效鉀、pH 變化明顯，種植1～20 a的露天菜地除速效鉀外其余土壤屬性均低于設(shè)施菜地。設(shè)施菜地種植時間越短土壤屬性指標值越高，而露天菜地種植時間越短全氮、速效鉀、速效磷的指標值越高，其他屬性指標值規(guī)律不明顯，露天菜地pH 在種植時間為20～30 a 時，指標值反而較高。從標準差可知露天菜地土壤屬性中，有機質(zhì)、全氮、速效鉀離散程度隨種植時間的增加而降低，pH 離散程度隨種植時間的增加而提高，有效磷離散程度規(guī)律不明顯；設(shè)施菜地土壤屬性離散程度隨種植年限的增加而降低。

表1 原始采樣數(shù)據(jù)分析Table 1 Analysis of raw sampling data

表2 不同城鎮(zhèn)距離土壤屬性指標平均值Table 2 Average value of soil properties in different town distance

2.1.2 菜地核密度空間分布

由圖3 可知，兩種菜地都具有空間集聚特征，露天菜地集中區(qū)域分布較廣，主要在鄉(xiāng)鎮(zhèn)和城市周邊（圖3 左圖）。設(shè)施菜地比露天菜地更加分散，主要集中于肥東縣西南區(qū)域（圖3右圖）。

2.2 菜地土壤屬性空間格局分析

如表4所示，全局Moran′sI值表明各土壤屬性均有顯著空間自相關(guān)性（P＜0.01），全局Moran′sI值排序結(jié)果是速效鉀＞有機質(zhì)＞全氮＞pH＞有效磷。

二元局域性莫蘭指數(shù)（Bivariate local Moran′sI）的結(jié)果表明土壤屬性和兩種菜地均有顯著空間自相關(guān)性（P＜0.01），與土壤有機質(zhì)、全氮及有效磷均呈空間正相關(guān)，與速效鉀和pH 呈空間負相關(guān)（表5），除速效鉀外，其余露天菜地的土壤屬性相關(guān)性均大于設(shè)施菜地。

二元局域Moran′sI可以定位雙因素空間自相關(guān)性，可確定菜地分布對周邊區(qū)域土壤屬性的影響類型，具體分為4 種：（1）High（土壤屬性值）-High（菜地核密度值）區(qū)域，代表菜地發(fā)展集中度高，土壤屬性累積與菜地相關(guān)；（2）High（土壤屬性值）-Low（菜地核密度值）區(qū)域，代表菜地發(fā)展集中度較低，土壤屬性累積較為嚴重的區(qū)域，但與菜地種植無關(guān)；（3）Low（土壤屬性值）-High（菜地核密度值）區(qū)域，代表菜地發(fā)展集中度高，無養(yǎng)分累積的區(qū)域；（4）Low（土壤屬性值）-Low（菜地核密度值）區(qū)域，代表菜地產(chǎn)業(yè)發(fā)展集中度低，土壤屬性累積較低的區(qū)域。如圖4 和圖5 所示，兩種菜地對5 種土壤屬性中的空間分布影響相同，有機質(zhì)和全氮在土壤中均呈現(xiàn)高高（High-High）空間正相關(guān)，有效磷和速效鉀均呈現(xiàn)低高（Low-High）負相關(guān)，pH呈現(xiàn)高低（High-Low）負相關(guān)。

表3 不同種植年限蔬菜種植類型土壤屬性平均含量的統(tǒng)計Table 3 Statistics of average soil properties content of different vegetable planting types and different planting years

圖3 露天菜地與設(shè)施菜地核密度分布Figure 3 Kernel density distribution in open vegetable fields and greenhouse vegetable fields

表4 不同土壤屬性全局Moran′s I值Table 4 Global Moran′s I values of different soil nutrients

3 討論

3.1 土壤屬性特征分析

研究區(qū)土壤有效磷的變異系數(shù)最高，速效鉀、全氮、有機質(zhì)及pH 緊隨其后。土壤屬性空間變異主要與人為生產(chǎn)或自然環(huán)境相關(guān)。因農(nóng)戶施肥習(xí)慣不同且磷在土壤中難以遷移造成其空間變異最高，鉀因容易形成半徑較小的離子且遷移性最強而空間變異位居其次，因土壤具有巨大的酸堿緩沖能力，因此土壤pH 變異系數(shù)最低，而有機質(zhì)與全氮遷移能力在速效鉀與pH之間，這與前人研究結(jié)果相似[29]。

各土壤屬性與城鎮(zhèn)距離的關(guān)系變化，源于土地利用變化，肥東縣最早發(fā)展蔬菜產(chǎn)業(yè)的區(qū)域一般位于近郊，隨著城市擴張，蔬菜產(chǎn)業(yè)種植中心發(fā)生位移，有效磷含量和速效鉀含量在2～2.5 km 處出現(xiàn)最高值。

露天菜地相比設(shè)施菜地土壤酸化明顯，這與前人研究結(jié)果相似[30]。露天蔬菜種植歷史相對較長，部分土壤酸化程度高于設(shè)施菜地，pH 相對較低，但種植時間最長的露天菜地因距離城市較近，隨著城市擴建所用建筑材料的影響，pH 值反而較高。設(shè)施菜地種植時間越短土壤屬性越高，因近10 a內(nèi)新增菜地多為土地流轉(zhuǎn)后的農(nóng)業(yè)合作社，為追求經(jīng)濟效益，其施肥量遠超過露天菜地，二者相差4～10倍，各土壤屬性除速效鉀外均高于露天菜地。露天菜地種植時間越短全氮、速效鉀值越高，其他屬性規(guī)律不明顯，這與農(nóng)戶長期形成的相同施肥習(xí)慣有關(guān)，但露天菜地地表容易產(chǎn)生徑流，具有較高的面源污染風(fēng)險。研究區(qū)露天菜地超20 a 種植區(qū)域集中在龍?zhí)拎l(xiāng)周圍，近10～20 a 的遍布肥東縣各鄉(xiāng)鎮(zhèn)周圍，由于設(shè)施種植是露天種植利潤的3～4 倍[31]，造成近10 a 肥東縣設(shè)施菜地興起，露天菜地增長較慢。10 a 以上設(shè)施菜地主要集中在牌坊回族滿族鄉(xiāng)內(nèi)，是因為將設(shè)施菜地作為扶貧重點項目發(fā)展，而其他鄉(xiāng)鎮(zhèn)發(fā)展緩慢。近10 a設(shè)施菜地增長迅速，遍布研究區(qū)所有鄉(xiāng)鎮(zhèn)，但集中于其西南部，這與研究區(qū)南部大開發(fā)后人口集中有關(guān)，但到達六家畈鎮(zhèn)后，因距巢湖太近，受到農(nóng)業(yè)環(huán)保限制而戛然而止。

3.2 空間格局分析

局部Moran′s I 指數(shù)對空間離群值敏感，可以對土壤中氮、磷等元素劃定區(qū)域，可為污染源的識別與控制提供依據(jù)[24]。一些研究表明該方法可得到各土壤養(yǎng)分元素含量分布的“高高”、“低低”聚集區(qū)和“低高”、“高低”區(qū)域的具體位置，可以分別從鄉(xiāng)鎮(zhèn)尺度[32]、行政村尺度[33]對土壤屬性進行管理。

本文與木合塔爾·艾買提等[34]的研究結(jié)果不同，后者結(jié)果中土壤pH 和緩效鉀的含量呈空間正相關(guān)，全氮在空間上不具有相關(guān)性，因本文以1 km 的網(wǎng)格作為評價單元，而后者以行政村為評價單元，導(dǎo)致空間尺度差異，從而造成研究結(jié)果不同。本文中設(shè)施菜地和露天菜地種植模式不同，對土壤屬性空間分布的影響具有一定差異性，但種植模式差異無法解釋土壤屬性所有的空間變異性。其他學(xué)者的研究也說明土壤屬性空間分布除受到地形、土壤類型等自然因素的影響之外，還受灌溉能力和施肥方式等人為因素的影響，且人為因素越來越強，越來越復(fù)雜[35]。

表5 露天菜地不同土壤屬性二元局域Moran′s I值Table 5 Bivariate local Moran′s I values

圖4 露天菜地土壤屬性二元局域Moran′s I分布圖Figure 4 Bivariate local Moran′s I distribution of soil nutrients in open vegetable fields

4 結(jié)論

（1）露天蔬菜種植歷史相對較長，距今10～20 a種植面積增長較快，主要集中于各鄉(xiāng)鎮(zhèn)周邊，種植時間越短其全氮、速效鉀的平均值越高，種植時間最長區(qū)域的pH受到城市開發(fā)的影響反而最高。

（2）近10 a設(shè)施菜地因利潤較高而種植面積增長迅速，種植時間越短土壤屬性指標值越高，除速效鉀外其余4 種屬性值均高于露天菜地。距離城鎮(zhèn)越近，土壤有機質(zhì)、全氮和pH的指標值越高，隨著距離的增加呈下降趨勢。

圖5 設(shè)施菜地土壤屬性二元局域Moran′s I分布圖Figure 5 Bivariate local Moran′s I distribution of soil nutrients in greenhouse vegetable fields

（3）兩種菜地種植與土壤酸化、土壤屬性累積相關(guān)性顯著，與有機質(zhì)和全氮呈現(xiàn)高高空間正相關(guān)，與有效磷和速效鉀呈現(xiàn)低高空間負相關(guān)，與pH 呈高低空間負相關(guān)。

（4）通過Moran′sI空間分析方法，可實現(xiàn)菜地對周邊土壤屬性空間格局影響的精準表達，實現(xiàn)菜地種植對不同區(qū)域影響差異的空間可視化，可以為進一步分析土壤屬性擴散演化機制提供參考，實現(xiàn)菜地不同區(qū)域及尺度的土壤屬性分區(qū)管理，有助于定義更嚴格的農(nóng)業(yè)種植管理規(guī)范，實現(xiàn)基本農(nóng)田的長期管理和保護。