东北黑土区典型县域种植模式遥感识别与空间格局分析

杜国明，张 瑞，梁常安，胡明宇

（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种类型，其中以“无序种植-大豆连作-两年轮作”类型为主，涉及的行政村广泛分布于在县域西北-东南方向；“大豆连作-无序种植-两年轮作”和“无序种植-玉米连作-两年轮作”种植模式结构次之，两者占比相当，前者在县域东部零散分布，后者呈东北-西南向带状分布；“无序种植-两年轮作-大豆连作”和“无序种植-两年轮作-玉米连作”两类涉及的行政村零散分布于县域东部和东南部。