東北黑土區典型縣域種植模式遙感識別與空間格局分析

杜國明，張 瑞，梁常安，胡明宇

（1. 東北農業大學公共管理與法學院，哈爾濱 150030；2. 東北農業大學經濟管理學院，哈爾濱 150030）

0 引 言

東北黑土區是全國重要的大豆和玉米主產區，在保障國家糧食安全方面具有舉足輕重的作用。黑土地雖然土地肥沃、糧食產量較高，但近些年來也面臨著土壤侵蝕加劇、黑土層變薄、土壤養分流失等生態環境問題[1]。2020年，習近平總書記在吉林考察時提出“保護好黑土地這一‘耕地中的大熊貓’，留給子孫后代”。作物種植模式對水土保持及黑土資源可持續利用具有重要影響，已有研究表明，相比于連作，豆米輪作可以明顯改善黑土土壤理化性狀，如提高土壤有機質含量、增加微生物群落多樣性和減少病蟲害影響等[2-5]，進而促進作物生長發育，實現大豆和玉米增產[6-7]，是一種能夠均衡用地養地、協調經濟效益和生態效益的種植模式。但在2007－2016年間，由于國家玉米臨儲和各項支農惠農補貼政策的實施，在市場作用下，東北地區糧食生產結構呈現“玉米增、大豆減”的時空演變基調[8]，傳統的“大豆-玉米”輪作種植模式受到嚴重破壞[9]。為促進農業可持續發展，國家于2016年出臺《探索實行耕地輪作休耕制度試點方案》，方案提出在東北冷涼區開展輪作試點，并建立了相應的輪作補助發放制度。因此，開展東北黑土區作物種植模式研究，對于探明區域輪作水平、優化輪作模式以及加強黑土地保護，推進“藏糧于地”戰略措施和質量興農戰略實施具有重要意義。

種植模式是指某一地區或某一耕地地塊在不同復種模式下作物的種植順序，主要包括連作模式和輪作模式，其中連作是指在同一地塊連季或連年種植相同作物的種植方式，而輪作則是指在同一地塊上有順序地在季節間或年際間輪換種植不同作物的種植方式[10]。多熟種植地區的種植模式是季節間的，由于復種模式復雜，學者們多從高時間分辨率遙感數據中提取時間序列植被指數，并根據植被指數的變化規律判斷作物種植順序[11]，進而揭示種植模式的空間格局特征。如張霞等[12]利用全年23個時相的MODIS EVI時間序列數據結合作物物候信息制作了華北平原兩熟制地區的種植模式分布圖。顧曉鶴等[13]從時間跨度為1.5 a的9幅遙感影像中提取NDVI時間序列信息，并采用變化向量分析法識別了北京順義地區的耕地輪作模式。許青云等[14]則利用460個時相的影像建立時間跨度為10 a的NDVI 長時間序列數據集和農作物季相節律特征識別了陜西省各年度的作物種植模式。郭昱杉等[15]基于MODIS 8 d NDVI時序數據識別出黃河三角洲地區的主要種植模式為冬小麥-夏玉米。Nguyen等[16]通過重建16 a期的MODIS EVI曲線以識別種植模式，進而劃分出越南中部Dak Lak省的單季和雙季玉米種植區。不同于多熟地區，一年一熟地區的種植模式是年際間的，低復種指數下僅需1個時相的影像即可實現作物類型提取，而后通過疊加方式獲取種植模式是學者們普遍采用的思路。如劉佳等[17]通過疊加2個相鄰年度的作物空間分布數據識別了黑龍江省海倫市的年際間作物變化信息。于鳳榮等[18]以黑龍江省的友誼農場為研究區，在3 a期作物分類數據的基礎上統計了各年度的作物面積，使用重茬率、迎茬率、重茬相對變化率等數理統計指標構建了作物重迎茬面積監測方法。總體來看，多熟地區和一年一熟地區種植模式的識別方法是有較大差別的，黑土區作為典型的一年一熟地區，部分學者已對該地區的種植模式監測進行了初步探索，但局限于2個相鄰年度間作物變化信息的識別和重迎茬面積的監測。一方面，種植模式不能簡單地等同于作物變化信息，它是后者的進一步分類和綜合，而這個過程是以豐富的作物變化信息為基礎的；另一方面，種植模式格局是在自然地理因素和社會經濟因素共同作用下形成的，除面積信息外，其空間分布信息亦是不能忽略的。因此，適當增加種植模式監測研究的時間跨度并使用能夠保留種植模式空間位置信息的方法對于揭示黑土區種植模式的類型及其空間格局是有必要的。

本文提出一種基于地學信息圖譜方法的作物種植模式監測方法，以黑龍江省克東縣作為研究區域，利用6個時相的遙感影像實現6 a間的作物分布提取，疊加之后進行圖譜重構以識別種植模式，最后使用核密度估計法分析了各類種植模式的空間格局特征，以期為優化輪作方案、制定輪作制度提供科學依據。

1 研究區概況

克東縣地處小興安嶺和松嫩平原的過渡地帶，介于126°01′～126°39′E，47°42′～48°17′N之間，隸屬于黑龍江省齊齊哈爾市，東至北安市，南鄰拜泉縣，西與克山縣毗鄰，北與五大連池市相望，如圖1所示。

克東縣下轄4鎮3鄉1農場，共包括99個行政村、4個農場、3個林場（為便于描述，以下統稱為行政村）。縣域內地勢丘陵起伏，緩坡漫崗，海拔高度147～411 m，屬于中溫帶大陸性季風氣候，年平均氣溫為2.4 °C，年平均降水量為526.5 mm，降水主要集中出現在7－8月，占全年降水量的49.9%。全縣土地總面積為208 324.29 hm2，其中耕地面積為147 084.53 hm2，墾殖率達70.60%。主要土壤類型為黑土、草甸土、暗棕壤、沼澤土等，其中，黑土覆蓋面積最大，占耕地總面積的75.61%。克東縣位于東北典型黑土區腹地，土質肥沃，適宜農業耕作，農業經營主體以農戶為主，作物熟制為一年一熟，主要種植大豆和玉米。克東縣種植業較為發達，2016年，全縣糧食作物播種面積為12.16萬hm2，全年糧豆薯總產量352 327 t，種植業產值達到128 089萬元，占全縣生產總值的30%。

2 數據來源與研究方法

2.1 數據源與預處理

根據東北地區主要農作物物候特征和已有研究[19]，本文選擇克東縣8月下旬或9月上旬的2012年Landsat7 ETM+和2013－2017年Landsat8 OLI共計10幅遙感影像作為主要數據源（表1）。除此之外，還收集了克東縣第二次土地調查年度更新數據庫、克東縣行政區劃矢量數據以及縣域30 m分辨率的DEM數據等作為輔助數據。

表1 遙感影像信息統計Table 1 Remote sensing image information statistics

本文使用ENVI進行作物分類，數據處理流程如下：1）使用DEM對影像進行幾何精校正，并在進行輻射定標、大氣校正、圖像融合、圖像鑲嵌等一系列預處理后，使用耕地矢量數據建立掩膜對影像進行裁剪；2）選取5、4、3波段進行標準假彩色影像合成，根據各類作物的光譜特征，以目視判別方式分別建立大豆、玉米、水稻以及其他作物（雜糧、薯類、漢麻等經濟作物）的解譯標志；3）基于解譯標志，建立作物分類訓練樣本點集。其中，大豆和玉米種植面積較大，所以各選取80個樣本點，將其隨機分為50個訓練樣本、30個驗證樣本。水稻和其他作物種植面積較小，所以各選取50個樣本點，將其隨機分為30個訓練樣本、20個驗證樣本。為提高解譯的精度，通過計算不同農作物樣本數據之間的可分離度來確定訓練樣本和驗證樣本的選擇是否合理。以2017年8月31日的影像為例，克東縣大豆、玉米、水稻和其他作物的訓練樣本可分度見表2。樣本可分離度的取值范圍為[0，2]，若可分離度大于1.9，表明不同作物樣本之間的可分離性較好，樣本的選擇較為合理；若可分離度小于1.8，則表明不同作物樣本之間的可分離性較差，樣本需要重新選擇；若可分離度小于1則樣本可以進行合并。從結果來看（表2），除了大豆和水稻的可分離度小于1.9之外，其他的農作物樣本之間的可分離性均較高，說明選擇的農作物樣本較為合理；4）使用監督分類方法對2012－2017年克東縣各類作物進行分類。

表2 訓練樣本可分離度統計Table 2 Statistics of training samples reparability

2.2 研究方法

2.2.1 地學信息圖譜

地學信息圖譜分析方法由陳述彭院士在20世紀90年代提出，它綜合了景觀圖的簡潔性和數學模型的抽象性，是一種經典的地理時空復合分析方法[20]。圖譜具有要素“圖形”和“譜系”的雙重性質，“圖”可以直觀地表達要素的空間分布特征，而“譜”則可以反映要素的過程變化信息，圖譜結合即可綜合展示要素的時空演變規律。借鑒已有研究[21]，地學信息圖譜中要素柵格單元的計算公式為

式中C為研究期間要素類型變化的單元屬性值；A為前一階段的要素單元屬性值；B為后一階段要素單元屬性值。

本文以2012－2017年間的6期作物分類柵格數據為基礎，首先在ArcGIS中對其進行重分類處理，將表征大豆、玉米、其他作物和水稻的柵格單元屬性值分別編碼為1、2、3、4，而后采用ArcGIS空間分析模塊中的“地圖代數”功能進行空間疊加，即可得到克東縣2012－2017年間的作物種植模式柵格數據。具體計算公式為

式中G為表征研究期間內作物種植信息變化特征的圖譜單元編碼值；Y2012、Y2013、Y2014、Y2015、Y2016、Y2017分別為表征2012、2013、2014、2015、2016、2017共6個年度農作物種植信息的圖譜單元編碼值。

2.2.2 核密度估計法

核密度估計法是在統計樣本數據的基礎上，通過獲取數據的概率分布曲線來擬合已知數據點，從而實現由已知估計未知的密度函數，屬于非參數檢驗方法。該方法通過計算要素在其周圍鄰域中的密度，可以反映點群的空間分布熱點、密度及其趨勢等。本文通過計算各類種植模式地塊的核密度，以展示種植模式的空間集聚特征。計算公式為

式中f(x)為核密度值；n為已知點數量；h為帶寬；k為核函數；x為估計點坐標值；xi為樣本點坐標值。

在核密度估計中，不同的帶寬h對密度分布結果有顯著影響，隨著h的增大，空間上點密度的變化更為光滑且概化程度更高；h減小時，估計點密度變化會顯的凹凸不平[22]。本文采用核密度估計法測度克東縣各類種植模式地塊的空間集聚特征。首先使用ArcGIS中“柵格轉點”功能將表征種植模式的柵格轉換為點狀要素，而后以此點狀要素為基礎進行核密度值的測算。

3 結果與分析

3.1 作物分類精度評價

由于目前國內尚無針對東北地區種植模式的遙感提取結果作為參考，無法直接評價本文種植模式的提取精度。本文以混淆矩陣的形式使用驗證樣本對作物分類結果進行精度評價，評價指標包括總體分類精度、Kappa系數、用戶精度、制圖精度4項，以2017年為例（表3）。克東縣2017年作物分類總體精度為92%，Kappa系數為0.89，大豆、玉米、水稻和其他作物的用戶精度分別為96.30%、93.33%、90.91%、85.71%，這表明分類精度滿足后續分析要求。

表3 克東縣作物分類混淆矩陣Table 3 Confusion matrix of crop classification in Kedong County

3.2 克東縣2012－2017年作物空間布局

根據作物分類結果對克東縣2012－2017年各類作物的種植面積進行統計，如圖2所示。

由圖2可知，2012－2017年，克東縣大豆與玉米的種植面積之和占總面積的比例均高于94%，水稻和其他作物的種植面積占比均較低且變化相對穩定。從大豆和玉米種植面積的年度變化來看，2013－2015年，大豆種植面積由9.20×104hm2下降至5.75×104hm2，玉米種植面積則由4.64×104hm2增加到8.28×104hm2。從2016年開始，玉米的種植面積顯著下降，由2015年的8.28×104hm2下降至2017年的3.74×104hm2，而大豆種植面積持續升高，由2015年的5.75×104hm2上升到2017年的1.01×105hm2，達到6 a間的最高值。

由于相鄰年度的作物空間分布情況相似，本文對主要作物種植面積的變化轉折點和極值點（2013年、2015年、2017年）的空間分布情況進行可視化處理（圖3）。可見，大豆、玉米種植地塊在縣域廣泛分布，6 a間兩者的分布區域不斷變化且呈“互為進退”關系；水稻在縣域北部呈東西向沿河流谷地帶狀分布；其他作物則在縣域內零散分布。

3.3 克東縣種植模式類型

使用地學信息圖譜方法，將6 a期作物分類柵格數據進行疊加之后，獲取作物變化信息圖譜。經統計，作物變化信息圖譜單元有2 857種編碼，依據種植模式的概念以及編碼特征進行分類[10]，即圖譜重構，大體可以歸并為6類種植模式（表4）。其中，僅由大豆和玉米兩種作物單元疊加產生的編碼有64種，可以發現：

1）表征某地塊在6 a間至少連續4 a種植大豆的編碼有8種，將其統一命名為“大豆連作模式”，此類模式意味著地塊的大豆連作問題較為突出；

2）表征某地塊在6 a間至少連續4 a種植玉米的編碼有8種，將其統一命名為“玉米連作模式”，此類模式則表示地塊的玉米連作障礙現象較為明顯；

3）表征某地塊在6 a間無明顯周期性豆米輪作規律的編碼有27種，本文將其統一命名為“無序種植模式”，此類模式在一定程度上實現了作物的交替種植；

4）表征某地塊在6 a間以3 a為周期，進行2個周期的以玉米為基礎茬口、大豆為中軸作物的“玉米-玉米-大豆”或以大豆為基礎茬口、玉米為中軸作物的“大豆-大豆-玉米”的種植模式編碼有6種，將其統一命名為“三年輪作模式”，這是耕地輪作的另一種形式；

5）表征某地塊在6 a間以2 a為周期，至少進行2個周期的“大豆-玉米”的交替種植的編碼有15種，將其統一命名為“兩年輪作模式”，此類模式是東北地區傳統的豆米輪作模式。

除此之外，共有2 793種編碼中至少含有一個代碼3或4，這表征某地塊在6 a間至少有1 a種植水稻或其他作物，將其統一命名為“混合種植模式”，其中包含水稻連作、水旱輪作、雜糧輪作等眾多種植模式。

根據種植模式分類模型，將作物變化信息圖譜中表征大豆連作模式、玉米連作模式、無序種植模式、三年輪作模式、兩年輪作模式、混合種植模式等6類種植模式的地塊進行制圖（圖4）。由于混合種植模式編碼類型較為繁雜，且面積占比（21.48%）相對較低，為簡化分析，本文以下內容對混合種植模式不做過多討論，僅選擇表4中種植模式二級類的前五項作為主要研究對象。克東縣各類種植模式的面積由大到小分別為無序種植模式、大豆連作模式、兩年輪作模式、玉米連作模式、三年輪作模式，其中，無序種植模式、大豆連作模式以及兩年輪作模式的面積占比之和為83.90%，是克東縣主要的種植模式。

表4 克東縣種植模式分類模型Table 4 Classification model of cropping patterns in Kedong County

3.4 克東縣種植模式的空間分布特征

本文采用核密度估計法測度克東縣各類種植模式地塊的空間分布。在帶寬的選擇上，以表征大豆連作模式的點狀要素為例，本文結合核密度分析工具中自動生成的1 600 m帶寬，分別設置1 000、2 000、3 000、4 000、5 000 m帶寬進行對比研究，據此比選最優帶寬。由于要體現縣域尺度上種植模式地塊的整體分布趨勢，因此選擇較大的帶寬，經過多次試驗，確定最優帶寬為4 km。在此帶寬下，既能較為清晰地辨別出各類種植模式地塊的密度中心，也能較完整地覆蓋各類種植模式地塊的分布區域（圖5）。

大豆連作模式地塊主要分布于縣域西部和北部，包括寶泉鎮、金城鄉以及千豐鎮；玉米連作模式地塊在縣域內呈現東北-西南走向的帶狀分布特征，其中以玉崗鎮最為明顯；無序種植模式地塊在縣域中部、東部和南部呈現集中連片分布的形態，在西部的金城鄉和北部的寶泉鎮亦有一定集聚分布；三年輪作模式地塊和兩年輪作模式地塊在縣域內均呈現“局部集聚、全域分散”的分布特征，同時，兩者在趙光農場北部的分布趨勢較為一致。將各類種植模式的分布區域進行對比，可以發現，無序種植模式地塊與大豆連作模式地塊存在一定的空間互補關系，與其他3類種植模式地塊的分布區域卻多有重疊，加之無序種植模式為縣域內面積占比最大的種植模式。因此，在縣域中部、東部和南部，無序種植模式可以視為克東縣種植模式的“基底”；對于面積占比為第二位的大豆連作模式地塊，與玉米連作模式地塊的主要分布區域呈現“互斥但不互補”的空間分布關系，與三年輪作模式地塊和兩年輪作模式地塊則亦有重疊。因此，在縣域西部和北部，大豆連作模式可以視為克東縣種植模式的“基底”。

3.5 克東縣種植模式的空間結構特征

為進一步分析各行政村的種植模式類型及比例，本文選取各行政村內面積占比為前3位的種植模式，并按其所占比例由高到低進行排列組合形成該村的種植模式結構，以展示行政村尺度種植模式的空間結構特征。通過計算各行政村前3位種植模式的面積之和占村域內所有種植模式面積的比例，得出克東縣所有行政村上述比例的總體平均值為88.29%，這證明此方法是適用的。為便于表達和制圖，將大豆連作模式、玉米連作模式、無序種植模式、三年輪作模式、兩年輪作模式分別賦予代碼“S、M、D、T、W”，最終將106個行政村歸并為10類種植模式結構（表5）。

表5 克東縣種植模式結構分區Table 5 Zoning of the cropping pattern structure in Kedong County

為便于進一步分析，選取行政村數量占比為前5位的種植模式結構（代碼為D-S-W、S-D-W、D-M-W、D-W-S、D-W-M）作為縣域的主要種植模式結構類型，并根據相似性原則將數量較少的S-D-T歸入S-D-W，D-S-M和D-S-T歸入D-S-W中，D-M-S和D-M-T歸入D-M-W中，可視化結果見圖6。由圖6可知，“無序種植-大豆連作-兩年輪作”種植模式結構（代碼為D-S-W）呈西北-東南向在縣域內連片分布，包含37個行政村，是縣域內占比最大的種植模式結構類型。該區域年際間的作物種植隨意性較高，且大豆連作障礙現象相對明顯，輪作以兩年輪作模式為主但在該區不占優勢，這與克東縣整體種植模式結構是一致的；“大豆連作-無序種植-兩年輪作”種植模式結構（代碼為S-D-W）零散分布于縣域西部和北部，在西南部則呈集聚分布特征，包含22個行政村。總體來看，該地區大豆連作問題比較突出，農戶種植行為亦缺乏引導，輪作政策實施效果較差；“無序種植-玉米連作-兩年輪作”種植模式結構（代碼為D-M-W）在縣域內呈東北-西南向帶狀分布，包含21個行政村。該區域作物種植隨意性亦較高，玉米連作障礙現象相對明顯，兩年輪作模式有一定分布但不占優勢；“無序種植-兩年輪作-大豆連作”種植模式結構（代碼為D-W-S）在縣域中部呈南北向離散分布，包含16個行政村。該區域的作物種植缺乏引導，但兩年輪作模式的推廣效果相對較好，大豆連作現象占較小比例；“無序種植-兩年輪作-玉米連作”種植模式結構（代碼為D-W-M）在縣域中部、南部和東部零散分布，涉及10個行政村。該區域作物種植無序特征相對明顯，兩年輪作效果較為突出，但玉米連作障礙現象同樣存在。

4 討 論

農戶的作物種植行為受農產品價格政策、市場價格、輪作種植等因素影響較大[9]。2013－2015年，克東縣大豆種植面積下降，玉米種植面積上升，這可能是受到大豆進口價格持續走低的影響，加之玉米單產及均價均高于大豆，導致農戶大豆種植意愿較低。從2015年開始，國家為調整糧食生產結構，鼓勵農民種豆積極性，在東北三省下調臨儲玉米掛牌收購價格，并調減“鐮刀灣”地區的玉米種植面積，造成玉米的種植面積顯著縮減，大豆種植面積快速攀升，調控政策富有成效。

克東縣無序種植模式占比最大，這可能是由于農戶作為“理性經濟人”，其作物種植決策行為易受農產品經濟收益影響而不斷調整；由于該地農戶已具備成熟的大豆種植技術和相對完善的機械設施，種植玉米反而會導致收益下降，故大豆連作現象較為嚴重；兩年周期的豆米輪作是當地傳統的輪作方式，農戶出于保護耕地的目的，使得兩年輪作模式的面積占比相對較高。

與劉佳等[17]的研究相比，本文在獲取作物變化信息的基礎上進一步識別出了種植模式。與于鳳榮等[18]的研究相比，本文除了分析種植模式的面積，還對種植模式的空間分布特征進行了測度。地學信息圖譜作為一種經典的地學分析方法論，在土地利用變化、景觀格局演變、生態信息表達等方面具有廣泛的應用[23]，但是在農業領域的應用研究相對較少。同時，種植模式是一種典型的時空復合體[24]，這與地學信息圖譜可以綜合表達空間位置和變化過程信息的內涵是一致的。因此，使用圖譜方法進行種植模式的識別是可行的。

為實現土壤改良、牧業發展、生態保護、增加糧油供給等目標，農業部、中央農辦等于2016年發布《探索實行耕地輪作休耕制度試點方案》，提出在東北冷涼區推廣“一主四輔”輪作模式，但未明確各類模式的具體輪作區域。因此，立足東北地區不同區域的資源條件和生態特點，探索各地區適宜的輪作作物、輪作格局、輪作周期并構建多目標協同的輪作區劃應當是下一步研究的重點[3]。

作物空間分布格局包含種植結構、種植面積等信息，是區域種植模式形成的基礎[25]。反之，在區域理想的輪作模式框架下，如何在空間維度上合理安排作物布局，在時間維度上協調作物種植順序，進而推動區域種植業結構調整和優化東北地區專業作物帶，對于確保國家糧食安全意義重大。

對于一年一熟地區種植模式的后續研究應注意選取適宜的時間尺度。過小的時間尺度難以識別出周期性的輪作規律，如2 a的時間尺度僅能體現作物變化信息；同時，農戶的作物種植決策行為亦在不斷調整中，多周期的輪作是難以實現的，過大的時間尺度下，種植行為的無序特征將愈加明顯。總體來看，時間尺度在3～6 a之間選取是較為合適的。

5 結 論

本文利用Landsat 8 OLI遙感數據提取克東縣作物類型，并結合地學信息圖譜方法識別種植模式，進而分析了種植模式的類型、面積、空間集聚性和空間結構等空間格局特征。主要結論如下：

1）2012－2017年，克東縣大豆和玉米種植總面積占比均高于94%，且兩者的種植面積在6 a間分別呈現“先減后增”和“先增后減”的變化趨勢。水稻和其他作物的種植面積占比較低且變化相對穩定。

2）克東縣各類作物種植模式按面積由大到小排序為無序種植模式、大豆連作模式、兩年輪作模式、玉米連作模式、三年輪作模式。其中，前3位種植模式面積占比之和為83.90%，是克東縣主要的種植模式。

3）大豆連作模式在縣域西部和北部占有明顯優勢，無序種植模式在縣域中部、東部和南部則是主要的種植模式。玉米連作模式呈東北-西南向帶狀分布，三年輪作模式和兩年輪作模式呈“局部集聚、全局分散”的分布形態。

4）行政村尺度的種植模式結構大體可以分為5種類型，其中以“無序種植-大豆連作-兩年輪作”類型為主，涉及的行政村廣泛分布于在縣域西北-東南方向；“大豆連作-無序種植-兩年輪作”和“無序種植-玉米連作-兩年輪作”種植模式結構次之，兩者占比相當，前者在縣域東部零散分布，后者呈東北-西南向帶狀分布；“無序種植-兩年輪作-大豆連作”和“無序種植-兩年輪作-玉米連作”兩類涉及的行政村零散分布于縣域東部和東南部。