王強(qiáng)，鄭夢(mèng)蕾，葉治山，楊善蓮，馬友華

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

1 材料和方法

1.1 研究區(qū)概況

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

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

1.2.1 蔬菜種植類型獲取

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

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

1.2.2 土壤采樣與分析

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

圖2 采樣點(diǎn)分布圖Figure 2 Distribution of sampling points

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

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

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

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

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

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

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

2 結(jié)果與分析

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

2.1.1 土壤屬性特征分析

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

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

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

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

表2 不同城鎮(zhèn)距離土壤屬性指標(biāo)平均值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é)果是速效鉀＞有機(jī)質(zhì)＞全氮＞pH＞有效磷。

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

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

表3 不同種植年限蔬菜種植類型土壤屬性平均含量的統(tǒng)計(jì)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ù)最高，速效鉀、全氮、有機(jī)質(zhì)及pH 緊隨其后。土壤屬性空間變異主要與人為生產(chǎn)或自然環(huán)境相關(guān)。因農(nóng)戶施肥習(xí)慣不同且磷在土壤中難以遷移造成其空間變異最高，鉀因容易形成半徑較小的離子且遷移性最強(qiáng)而空間變異位居其次，因土壤具有巨大的酸堿緩沖能力，因此土壤pH 變異系數(shù)最低，而有機(jī)質(zhì)與全氮遷移能力在速效鉀與pH之間，這與前人研究結(jié)果相似[29]。

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

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

3.2 空間格局分析

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

本文與木合塔爾·艾買提等[34]的研究結(jié)果不同，后者結(jié)果中土壤pH 和緩效鉀的含量呈空間正相關(guān)，全氮在空間上不具有相關(guān)性，因本文以1 km 的網(wǎng)格作為評(píng)價(jià)單元，而后者以行政村為評(píng)價(jià)單元，導(dǎo)致空間尺度差異，從而造成研究結(jié)果不同。本文中設(shè)施菜地和露天菜地種植模式不同，對(duì)土壤屬性空間分布的影響具有一定差異性，但種植模式差異無法解釋土壤屬性所有的空間變異性。其他學(xué)者的研究也說明土壤屬性空間分布除受到地形、土壤類型等自然因素的影響之外，還受灌溉能力和施肥方式等人為因素的影響，且人為因素越來越強(qiáng)，越來越復(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）露天蔬菜種植歷史相對(duì)較長，距今10～20 a種植面積增長較快，主要集中于各鄉(xiāng)鎮(zhèn)周邊，種植時(shí)間越短其全氮、速效鉀的平均值越高，種植時(shí)間最長區(qū)域的pH受到城市開發(fā)的影響反而最高。

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

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

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

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