基于像元尺度耕地質量局部空間自相關的基本農田劃定

劉彥文 劉成武 何宗宜 周 霞

(1.武漢大學資源與環境科學學院， 武漢 430079; 2.湖北科技學院資源環境科學與工程學院， 咸寧 437100; 3.中南民族大學公共管理學院， 武漢 430074)

0 引言

農田為人類糧食安全提供了保證，劃定永久性基本農田將保障中國農業發展的生產基線[1]，我國農田資源的大量流失引起了各級政府部門的關注[2]，近年來劃定了基本農田保護區、糧食生產功能區和重要農產品生產保護區，對優質農田進行永久性保護。基本農田劃定是一個復雜的過程[3]，雖2017年完成了全國1.03億hm2永久基本農田的劃定工作，且較以往劃定成果質量有明顯提升[4]，但由于優質耕地選擇標準較難統一，劃定成果質量仍有改進空間。為此，需探討相關劃定方法，為實現數量、質量、生態“三位一體”耕地保護新格局，制定后續兩區科學合理的劃定標準，具有重要意義。

目前研究對象以基本農田、高標準基本農田、永久基本農田等為主；研究內容主要集中在從耕地質量、土地利用、空間布局等出發的農田保護、劃定/分區、結果評價等方面；GIS是主要研究方法，也有層次分析法、局部空間自相關、LESA、TOPSIS、熵權法、理想點法等。在制定優質耕地入選為基本農田的標準時，基于耕地質量的評價體系占有較大比重。在評價指標體系中，耕地自然質量、區位條件最為常用[5-11]。有關學者考慮耕地空間布局與社會發展因素，引入空間形態、社會經濟等相應指標對指標體系進行完善[12-14]。耕地兼有生態服務功能，生態環境與耕地穩定性指標也被相關學者關注[15-17]。在研究尺度方面，多以市縣級為對象，從耕地圖斑角度出發進行分析[18]，從柵格網與像元粒度的分析則更加細致[19-20]。鑒于耕地在空間上分布存在一定的規律性，文獻[18,21-23]對基于局部空間自相關的方法進行了研究。

基本農田必定是耕地，但耕地不一定是基本農田，選擇優質耕地的標準除了常用的耕地自然條件、區位條件等之外，可加入反映耕地利用情況、作物長勢等情況的內因指標。本文以關聯反映耕地利用狀態與作物自身信息的NDVImax-min、NDVIsd作為完善耕地質量的綜合評價指標。以湖北省嘉魚縣為研究對象，將所有資料歸一化至30 m×30 m柵格像元尺度，在采用局部自相關分析像元耕地質量綜合指標值的基礎上，對基本農田進行劃定，以期為進一步完善耕地優選為基本農田的指標體系、提升縣級基本農田劃定綜合質量評價至像元級提供相關參考。

1 研究區與數據來源

1.1 研究區概況

嘉魚縣位于湖北省東南部、長江中游南岸(圖1)，地理坐標為113°39′～114°22′E，29°48′～30°19′N。縣域地處江漢沖積平原，屬副熱帶濕潤季風氣候，日照充足，雨量充沛，土地肥沃，全縣轄8個鄉鎮1個農場，人均土地0.28 hm2，人均耕地約0.09 hm2。地勢自西南丘崗向東北平原交迭過渡，總體海拔較低，大部分在19～50 m之間，全縣大體形成“一山三水四分田、兩分道路與莊園”的地貌格局。

圖1 研究區示意圖Fig.1 Location of study area

1.2 數據來源與處理

數據主要包括國土部門提供的土地利用現狀、土地利用規劃、基本農田和土地利用分等定級等數據庫資料；研究區域土地利用現狀、市縣鄉村四級土地利用意見等調研資料；從地理空間數據云(http:∥www.gscloud.cn/)下載的研究區自然年內15期Landsat遙感影像，購買研究區GF1影像，研究區社會經濟數據以及Google Earth影像數據；通過實地考察與調研進行資料分析判斷與綜合取舍，與已有相關數據、互聯網收集資料等進行交叉驗證，以檢驗數據來源的真實可靠性。將所有數據按照評價指標體系歸類，在ArcGIS中將每個指標數據處理為一層柵格數據[24]。

2 研究方法

2.1 耕地質量綜合評價

2.1.1評價單元與指標體系

縣級是基本農田劃定的落實單位，評價單元是劃定成果構成的基本單元，其尺度與劃定成果質量直接相關。在市縣級尺度上的評價單元包括耕地圖斑[18, 22]、不同尺度柵格單元[19,24]等，本文以30 m×30 m柵格為基本評價單元。

《基本農田劃定技術規程(TD/T1032—2011)》是基本農田劃定的基本依據，結合區域農用地分等成果、區域基本農田劃定細則等，現有指標體系基本覆蓋耕地自然質量、耕地立地條件、耕地利用、社會經濟和生態環境條件等基本準則。以農田作物自身為參照，可以將現有指標體系歸為外部因素，其反映了外界環境對劃定成果的影響。所有的外部因素最終反映在農田自身方面，植被NDVI時序數據記錄了不同時段作物自身的狀況，文獻[25-26]對NDVI時序數據與耕地利用狀況進行了相關研究，同一耕地空間單元如在自然年內NDVI發生了變化，說明對應的土地被實際利用過，NDVI時序數據的統計變化特征可以輔助完善基本農田劃定指標體系。本文選擇自然年內像元NDVI時序數據中的最大值與最小值之差NDVImax-min間接反映耕地利用與作物自身的差異，用NDVI數據序列的方差NDVIsd間接反映耕地利用頻率等的差異，并將其歸入耕地景觀生態條件，最后所建立的耕地質量綜合評價指標體系如表1所示。

表1 嘉魚縣耕地質量綜合評價指標體系Tab.1 Comprehensive quality evaluation system of cultivated land in Jiayu County

注：Pi=C表示將P指標對應變量i賦值為C，Ed指歐氏距離。

2.1.2評價指標賦值方法

對于距離連續型變量，首先運用ArcGIS中Euclidean Distance計算縣域內每個像元到對應目標的距離，然后視距離指標與優選基本農田準則的正負相關性，將所有值都歸一化在[0，100]以內。對于常規離散型指標，Q1、Q4賦值計算公式為

F+=100(Xi-Xmin)/(Xmax-Xmin)

(1)

F-=100(Xmax-Xi)/(Xmax-Xmin)

(2)

Q2、P1、P2參考《農用地質量分等規程》及已有文獻如表1賦值，pH值的賦值計算公式為[27]

(3)

其中a1～a4分別為5.0、6.5、7.5和8.5，Abs(·)為取絕對值函數。

像元連片度計算，針對每個像元，運用Qeen空間關系鄰接規則，判斷其周邊8個相鄰方向上是否有基本農田像元，每發現一個相鄰的基本農田像元，其像元連片度增加1。單個像元的連片度最小為0、最大為8，數值越大，像元連片度越好。優先保留原有基本農田中的高等級耕地、集中連片耕地，要求劃定后基本農田集中連片程度有所提高[28]。以耕地圖斑為基本單元進行研究劃定時著重關注圖斑的連片性，由于圖斑的平均尺度一般相對像元(約900 m2)而言較大，劃定精度不及像元尺度，所以本文以像元為基礎，逐像元計算其連片程度，最后再運用式(1)將像元連片程度進行歸一化。

2.1.3網絡層次分析法確定指標權重

耕地質量綜合評價指標部分相互獨立、部分存在關聯，網絡層次分析法(ANP)不但具有層次分析法(AHP)遞階式層次結構，而且還具有內部依賴性和反饋性的層次結構[29]。采用ANP與專家打分法相結合的方式確定指標權重。依據ANP各指標權重求解過程，首先對指標關聯情況進行分析，并構建判斷矩陣；其次，將專家評分結果輸入SuperDecision軟件進行一致性檢驗，對不滿足一致性檢驗的結果進行調整并反饋給專家再次征求意見；最后，輸入評分結果、求解超級矩陣計算指標權重(表1)。

2.1.4像元耕地質量綜合分值計算

像元是本文的基本評價單元，采用多因素綜合評價法對評價單元進行綜合質量評分，分值越高耕地質量越好[18]，像元綜合分值計算公式為

(4)

式中Zi——第i個像元的耕地質量綜合得分

wj——第j個指標的權重

Zij——第i個像元第j個指標的初始值

gstan(·)——標準化函數，保證像元在每個指標的得分都在[0,100]之間

2.2 局部空間自相關分析

近年來，很多學者基于耕地質量指數采用局部空間自相關方法對耕地合理利用進行了相關研究[21-23,30-31]。局部空間自相關指標(LISA)可以揭示空間參考單元與其鄰近單元屬性特征值之間的相似性或相關性，運用LISA指標中最常用的局部莫蘭指數(Local Moran’sI)和Moran散點圖對耕地質量在空間上的聚集、異質或隨機的分布特征進行研究，局部莫蘭指數計算式為

(5)

(6)

(7)

式中Ii——單元i的局部莫蘭指數

wij——單元i與j之間的空間權重

n——與單元i相鄰接的空間單元數量

xi、xj——單元i、j的耕地質量指數

σ2——方差

空間權重矩陣是進行局部莫蘭指數計算的關鍵，不同鄰接規則得到的空間權重矩陣不同，Rook、Queen是最常用的鄰接規則。耕地的連片性在中心像元周邊8個方向上均有體現，本文選擇一階Queen鄰接規則構建空間權重矩陣。局部莫蘭指數高值表明有相似變量值的面積單元在空間集聚，低值表明不相似變量值的面積單元在空間集聚，Moran散點圖將空間單元的集聚形式分4個象限進行可視化。第1、3象限分別為HH與LL正相關型，各自表示高值與高值、低值與低值集聚，第2、4象限分別為LH、HL負相關型，空間異質性主要區別是高值包含低值異常和低值包含高值異常。

2.3 耕地利用與像元NDVI時序數據關聯分析

NDVI時序數據在反映植被生長季節性和年內特征方面具有良好的可信性[32]，文獻[25-26]利用其對耕地、休耕地進行了相關研究，并取得了良好效果。同一像元由于自然年內不同種值季數的差異(如一年一季、一年兩季、一年多季等)會在NDVI時序數據中有不同取值，NDVI時序曲線峰值個數及其方差、最大最小值及其極差等都與耕地利用狀況相關聯。理論上曲線形態與種植收割相對應，裸地從作物種植開始到其生長最旺盛，從成熟到收獲，對應曲線的不同取值，曲線方差越小說明耕地種植利用越平穩。極值和極差可反映一個種植收割期內不同作物間的橫向生長狀況。利用NDVImax-min、NDVIsd指標關聯分析耕地利用狀況，其計算公式為

NDVImax-min=max(NDVIseries)-min(NDVIseries)

(8)

NDVIsd=SD(NDVIseries)

(9)

式中NDVIseries——NDVI時序數據序列

其中，max(·)、min(·)、SD(·)分別為針對NDVI時序序列取最大值、最小值和方差函數，數據越密集效果越好。

3 結果與分析

3.1 耕地綜合質量指標處理結果與分析

本文像元面積約900 m2，與耕地圖斑相比其單元面積較小且均一，同時與研究區人均耕地866 m2也最為接近。基本農田的劃定要求上圖入庫、落實到戶，像元尺度可在更細粒度上落實基本農田的保護范圍與保護責任。將評價耕地質量所需的16個指標在ArcGIS中各自處理成一幅30 m×30 m的柵格數據，統一采用WGS-1984-UTM-Zone-49N投影坐標，并采用Snap Raster保證不同指標層單元匹配對齊。全縣范圍內所有指標進行歸一化處理，范圍在[0,100]之間，結果如圖2所示，運用式(4)可得各像元耕地質量綜合指數。

圖2 像元尺度耕地質量評價指標柵格圖Fig.2 Raster maps of farmland quality evaluation indexes at pixel scale

由圖2可知，數值離散型指標(如Q1、Q2)由于原始評價基于耕地圖斑而賦值，其局部區域內取值較為集中，除了邊界更替處對像元區分有貢獻之外，同質區域內部對像元區分無貢獻。對于數值連續型指標(如L1、E4)，以像元為單元計算距離，可以在圖斑內部進一步細化距離權重對劃定結果的影響。相比以圖斑為單元計算圖斑質心與目標之間的距離、圖斑內部處處同值而言，基于像元的分析可以使劃定結果更趨合理，本文有7個指標與距離相關，占總指標的43.75%。

3.2 局部空間自相關結果與分析

采用GeoDa與ArcGIS軟件進行局部空間自相關分析與制圖，首先將像元轉為矢量點，并將像元耕地綜合質量增加到點屬性，其次將GeoDa計算結果保存到點屬性中，最后將矢量點帶屬性轉為與已有指標同樣的柵格進行可視化分析。

全縣耕地質量局部空間自相關指數為0.864 5，由于耕地像元與非耕地像元的耕地綜合質量評價得分呈兩級分化趨勢，所以Moran散點圖(圖3a)呈兩個團聚狀，綜合得分對耕地與非耕地具有良好的區分度。就耕地像元而言，其耕地質量局部空間自相關指數為0.991 6，其空間聚集程度較全縣更加明顯(圖3b)。

圖3 像元耕地質量空間關聯分析Fig.3 Spatial correlation analysis of farmland quality in pixel scale

由圖3c、3d可知，像元耕地質量空間關聯集聚效果明顯，在空間上，全縣范圍的HH型主要分布在東北平原區，包括牌洲灣鎮、潘家灣鎮、渡普鎮，以及新街鎮的大部分范圍，該區交通便利、自然條件優越，是典型的水稻蔬菜基地[33]，LL型在長江、斧頭湖等水域區最為明顯，縣域西南部以丘陵崗地為主，長江沿線南門湖村、官洲村等HH相對明顯外，主體區內其余類型交替相間，集聚規律不明顯。由圖3d可知，正相關類型與耕地質量關聯情況更為明顯，尤其體現在交通便利、地勢平坦、土質肥沃的省道S102兩側，但牌洲灣鎮北部及其與潘家灣鎮連接處相關性不明顯，可能受洪澇災害及近年來新起的畈湖工業園區的影響等所致。中部空白處是縣城所在地，其西側南門湖村是唯一保留的大面積耕地區，其HH型尤為明顯。負相關類型集中體現在官橋鎮、高鐵嶺鎮以及陸溪鎮的大部分耕地中，該區較東北部耕地分布較為零星、地勢錯落起伏、林地比重較大，像元耕地綜合質量較低，LL型聚集與此相對應。

由表2可知，在全縣范圍內HH與LL正相關類型占72.5%，負相關類型共占1%，且其中90%主要分布在HH型周邊、并與其形成了良好的連片分布。在耕地范圍內，正相關類型占59.3%，說明耕地質量在空間上集聚特征明顯。無顯著相關的區域占比40.7%，P1、E3與Q1、Q3等因子在此處高低值交替相間，空間關聯錯綜復雜。對耕地范圍內全縣集聚情況分析可知，HH型與耕地重合度為81.3%，說明即使在耕地范圍未知的情況下，若優先利用HH型進行基本農田劃定，也可達到較好的效果，這對于利用遙感反演信息客觀進行基本劃定提供了良好的依據。

表2 局部空間自相關結果分析Tab.2 Analysis of local spatial autocorrelation results

注：P<0.05，在95%置信度時統計結果；NS型表示非顯著型。

3.3 基于像元局部相關性的基本農田劃定結果與分析

局部空間自相關值Moran’sI反映了中心像元與周邊像元的相關程度，其在[-1，1]之間，越靠近1表示中心像元屬性與周邊像元同質屬性正相關性越強，反之，越靠近-1表示其負相關越強，等于0時表示規律性較弱，呈空間隨機分布狀態。依據其相關關系，相關學者按HH型、LL型、LH型和HL型從宏觀層面定性對耕地利用進行了人為劃分區塊并分析建設時序的研究，如呈空間擴散效應的HH型為高質量耕地圖斑被同樣高質量圖斑所包圍，所以應優先劃定、區域內禁止非農建設。本文從像元微觀層面出發，利用相關屬性值由大到小循環迭代依次選入單個像元，并將劃定結果與空間集聚情況進行聯合分析，具體實現過程利用Python完成。

全縣劃定任務為304 843個像元單位，由于在全縣范圍內土地利用類型復雜、質量評定指標多樣，局部空間自相關莫蘭指數理論上存在非耕地高于耕地的情況(如部分水域)，所以，選擇以下兩種方法進行劃定，方法1是在全縣范圍內按前述計算的像元耕地綜合質量得分，由高到低依次選入像元確定劃定區域；方法2是在耕地范圍內按莫蘭指數由大到小順序依次選入像元確定劃定區域，結果如圖4所示，對比分析如表3所示。由圖4a可知，不管是利用全縣范圍內耕地質量還是耕地范圍內莫蘭指數，其劃定結果都保持了高度的一致性，二者重合像元數為259 040個，占劃定任務的85.0%。結合耕地質量評價數據與Google高分影像疊置分析，基于耕地質量的劃定，其優勢主要在于劃入的像元是現有指標下耕地區內質量最優的區域，且與二者共有部分的連片性保持較好，基本在其周圍分布，但其對丘陵區隴形地帶等效果欠佳。基于莫蘭指數的劃定結果對縣域西南部丘陵崗地區存在的條形、隴形相對碎片耕地的區域劃定效果占優，但整體而言連片性不及前者。

由圖4b與表3可知，兩種方法都優先將空間上呈HH型聚集的99.97%像元劃入了保護范圍，這是劃定成果的重點保護區域，區域內應該禁止非農建設[21]。HH型中未劃入的29個像元雖然耕地質量較好，但大部分呈獨立像元形式存在，小部分僅一邊與HH相連，極大地影響了片塊景觀布局形態，在相關項目占用耕地時，可以考慮優先利用。兩種方法的差異主要在于對LL型和NS型的處理上，基于耕地綜合質量的劃定優先將質量得分較高但指數得分低、在空間上呈NS型聚集的部分劃入，其主要分布在牌洲灣鎮北部、畈湖工業園周邊、渡普鎮與斧頭湖連接處等。基于莫蘭指數的劃定則優先考慮了指數得分，與前者差異在于質量得分相對前者較低，在空間上呈LL型。此部分主要分布在陸溪鎮西部、高鐵嶺鎮東側，后者明顯可見丘陵崗地內隴形耕地輪廓特征。由于基本農田的劃定是一個復雜多層次的過程[34]，二者差異部分也體現出了各自的優勢，前者在片塊連接性方面相對較好，后者在丘陵崗地區劃定較為接近現實，二者空間上的直觀分布，為最終決策提供了良好參考。

圖4 劃定結果Fig.4 Graph of demarcation results

類型LISA莫蘭指數/耕地質量重合部分莫蘭指數獨立部分耕地質量獨立部分像元像元比例/%像元比例/%像元比例/%比例/%HH型11515311512432.90029032.9LL型926454695813.44566113.026026.4LH型7800007800HL型100000010000NS型1426709695827.714204557013.040.7小計35064625904074.04580313.04580313.0100

4 結論

(1)在縣級空間尺度內耕地綜合質量涉及距離類指標時，一般基于像元的評價較以耕地圖斑為單元的評價更具良好的區分性，劃入的基本農田理論上更趨合理，對人均、戶均基本農田保護責任的落實也更加精確。但基于像元的相關分析涉及數據量較大，對運算資源的要求也相對較高。

(2)湖北省嘉魚縣耕地綜合質量在空間上表現出較強的空間相關性，全域和耕地范圍的局部空間自相關指數分別為0.864 5和0.991 6。基于基本農田包含于耕地的先驗知識，運用后者的劃定更趨合理，其HH正相關類型主要集中在潘家灣鎮、渡普鎮等東北平原主導區，LL負相關類型主要集中在官橋鎮、高鐵嶺鎮等西南部丘陵地帶，且分布較為零星。

(3)基于耕地質量與莫蘭指數的劃定結果在空間上保持了高度的一致性，兩種方法都優先將空間上呈HH型聚集的99.97%的像元劃入了保護區，前者在綜合質量方面占優，但對丘陵區隴形地帶等劃定效果不及后者。基于像元的分析對于嘉魚縣的劃定更為合理，NDVI等指標間接關聯了土地利用狀態與農作物自身信息，對耕地質量綜合評價體系的完善有益。