基于NDVI时间序列影像的农业种植结构提取——以河北省邯郸市广平县为例

Extraction of Agricultural Planting Structure Based on NDVI Time Series Image—A Case Study of Guangping County, Handan City, Hebei Province

摘要: 农业种植结构提取是基础性环节，针对原来以人工传统调查统计的农作物信息获取方法主观性较强且误差巨大难以满足信息时代背景下对农业的管理。近年来，水资源匮乏不断制约着灌溉农业发展，不同作物在全生育期所需的灌溉量存在巨大差异，因此快速准确的了解种植结构对于节水农业有巨大效益。本文以广平县为研究区，以哨兵2号遥感影像作为基础数据，以人工采集和目视解译的样本点为标准，采用随机森林算法对广平县的种植结构进行分类。结果表明：分类结果的总体精度为89.95%，Kappa系数为0.88，总体精度和Kappa系数均高于0.85，说明各个农作物分类结果较为精准。经统计，广平县内冬小麦–夏玉米、棉花、玉米、红薯、谷子和花生的种植面积分别为148.23 km²、9.34 km²、18.25 km²、2.31 km²、0.77 km²和3.98 km²。

Abstract: Extraction of agricultural planting structure is a fundamental process. The traditional manual survey and statistical methods for obtaining crop information were subjective and prone to sig-nificant errors, which could not meet the management requirements of agriculture in the context of the information era. In recent years, the scarcity of water resources has constantly hindered the development of irrigation agriculture. Different crops have significant differences in the amount of irrigation required throughout their growth period. Therefore, a rapid and accurate understanding of the planting structure is of great benefit to water-saving agriculture. In this study, Guangping County was selected as the research area, and Sentinel-2 satellite remote sensing images were used as the basic data. With the standard of manually collected and visually interpreted sample points, the random forest algorithm was employed to classify the planting structure in Guangping County. The results showed that the overall accuracy of the classification results was 89.95%, with a Kappa coefficient of 0.88. Both the overall accuracy and Kappa coefficient were higher than 0.85, indicating that the classification results of various crops were accurate. According to statistics, the planting areas of winter wheat-summer maize, cotton, maize, sweet potato, millet and peanut in Guangping County were 148.23 km², 9.34 km², 18.25 km², 2.31 km², 0.77 km² and 3.98 km² respectively.

1. 引言

农业种植结构是计算农业产值、需水量、粮食产量的基础数据，为了能够保障农业的稳定发展，国家必须制订适合的农业政策和经济发展规划，这就需要准确及时的了解农业的种植面积 [1] 。由于过去科学技术发展的限制，农作物种植信息是通过人力统计的方式实现，这种方式需要层层上报，往往耗时耗力，而且在准确性时效性方面存在很大缺陷 [2] 。随着遥感技术的兴起，为农作物种植面积及其空间分布提取提供了新思路，Jakubauskas [3] 于2002年利用光谱信息的变化区分玉米大豆等作物，但分类结果欠佳。

之后，遥感解译工作逐步趋向成熟，同分辨率条件下，由单时相提取逐步变为多时相提取，分类精度有所提高。Wardlow B D [4] 于2008年利用NDVI和EVI两种时间序列识别大面积农田，并比较其差异得到最佳区分结果图；闫慧敏 [5] 等于2005年利用多时相遥感影像，对中国农田种植信息进行分类；Julien [6] 等于2011年采用最大似然法提取西班牙部分区域种植结构的变化，分类效果较为理想。同年，国内的赵丽花 [7] 、潘耀忠 [8] 利用EVI时间序列进行分类识别，得到较好结果。张健康 [9] 于2012年运用多时相的遥感影像数据对黑龙港地区的棉花、玉米、小麦和蔬菜进行农业种植结构的提取，结果显示其总体分类精度达到了91.3%，分类精度较高，能够较好的反映作物的分布情况，为该地区主要的种植结构调整提供科学依据。刘焕军 [10] 于2015年建立棉花的光谱信息变化时间序列图，实现对其产量的预测。

随着以美国国家航天局的Landsat系列等卫星发射和完善，其后随着中高分辨率卫星不断发射，尤其是以高分6号，哨兵2号等卫星为代表的开源卫星数据，使得分类精度再次突破。熊元康 [11] 等于2018年对天山北坡经济带的棉花、玉米、小麦利用随机森林分类器进行分类，总体精度为92.19%，Kappa系数为0.88；刘嘉文 [12] 于2020年以县域为尺度使用哨兵2号影像数据对桦南县主要农作物玉米、大豆水稻进行空间种植结构提取研究，总体分类精度95.22%，Kappa系数为0.89；2022年，田鑫 [13] 以甘肃张掖灌区为研究区，将基于时序NDVI特征和光谱特征提取的玉米面积相比，发现基于时序NDVI特征的提取的结果优于基于光谱特征提取的结果；2023年，牛昱杰 [14] 采用哨兵2号遥感影像，并基于多时间序列的归一化植被指数构建模型，获取了2018~2021年甘肃省甘谷县冬小麦的分布情况，提取的总体精度均在90%以上，为当地冬小麦遥感种植监测提供借鉴和参考。

河北省处于华北平原西北部是我国重要的农业基地，但人均水资源量仅为全国平均水平的11%，其中农业用水占比超过70%，且地下水占比的份额较大，是典型的资源型缺水省份 [15] 。从上个世纪70年代到现在，年均地下水超采量约为59.7亿m³，尤其是中东部的平原区已经成为超采最严重地区，近三十年超采量达1500亿m³，随着地下水位的持续下降，目前河北省地下水超采问题十分严重，爆发了一系列连带问题，上世纪六十年代仅邯郸市就有13条地裂发生，对当地的居民生活产生极大危害，并制约着社会经济的可持续发展 [16] 。邯郸市广平县作为河北省水资源超采区的农业县，摸清该县的种植结构无疑为水资源管理、粮食宏观调控和促进农业可持续发展提供有效途径。

2. 数据及方法

2.1. 研究区概况

本文选择广平县为研究区，该地区地处河北省南部，邯郸市东部，介于东经114˚84'至115˚19'，北纬36˚37'至36˚62'之间，全县境域长31.4 km，宽27.0 km，总面积达320 km²。县政府驻地广平镇，距省会石家庄198 km，距邯郸市中心城区45 km。全县地势起伏不大，自西南向东北逐渐走低，平均坡度1/3500，海拔高度在40~53.9 m之间，最大高差13.9 m。广平地理位置如图1所示。

Figure 1. Geographical location map of Guangping County

图1. 广平县地理位置图

2.2. 数据获取

本文选取哨兵2号(Sentinel-2)卫星作为原始数据，该卫星源于欧洲委员会和欧洲航天局共同倡议的哥白尼计划，由Sentinel-2A和2B两颗相同的卫星组成，其光学任务主要是对于陆地和沿海地区进行各种自然资源灾害的监测，涵盖可见光、近红外线与短波红外线的13个波段，空间分辨率为10~60 m，时间分辨率为5天。本文采用空间分辨率为10 m的B2、B3、B4、B8波段进行后续分析，共选取从2019年9月至2020年10月的27景L1C级无云影像数据。

2.3. 方法

2.3.1. 典型作物选取

根据当地调查的物候情况，以及种植面积的大小，选取冬小麦、夏玉米、棉花、花生、谷子、红薯共6种农作物作为本次种植结构提取的主要作物，其中由于冬小麦和夏玉米轮作情况较为普遍，故把本次遥感提取的农作物分为冬小麦–夏玉米、棉花、玉米、红薯、谷子、花生。

2.3.2. 图像预处理与NDVI计算

对27景原始遥感影像进行影像预处理操作。由于哨兵2号L1C级数据已正射校正和几何精校正，所以影像预处理主要包括辐射定标、大气矫正、影像裁剪三个部分。

辐射定标：将遥感影像传感器纪录的无量纲DN值转换成具有实际物理意义的表观辐射亮度的过程，其意义在于尽可能地除去由传感器自身造成的影响，减少遥感影像与地面物体辐射亮度之间的差异，使信息更加准确，为农作物种植结构提取、面积估算提供精准的数据支撑。计算公式如下 [17] ：

$Le=Gain\times DN+offset$ (1)

式中：Le为影像的辐射亮度值，Gain表示传感器定标系数的增益值，DN为影像的数字量化值，offset为定标系数的偏移值。

大气校正：太阳辐射与物体表面相互作用时，需要透过大气层，而大气中的气溶胶、二氧化碳、臭氧等会吸收部分反射光，同时伴随着少量大气散射现象，使得成像过程中农作物反射光谱信息出现了偏差，存在一定程度上的辐射变形，消除以上干扰的过程通常称为大气校正。

图像裁剪：图像裁剪是按照某一规定的边界进行裁剪，目的是去除研究区范围以外的多余部分，以免操作过程中带来的巨大工作量。

对预处理之后的27景影像使用归一化植被指数NDVI作为分类特征指数 [18] 。它由Rouse等提出，在许多研究中得到应用 [19] [20] 。该指数通过建立两波段间的数学运算关系，表达了红光波段与近红外波段即植物吸收强烈的波段与植被反射强烈的波段之间的关系，其计算公式如下 [2] ：

$NDVI=\frac{{\rho }_{NIR}-{\rho }_{Red}}{{\rho }_{NIR}+{\rho }_{Red}}$ (2)

式中： ${\rho }_{NIR}$ 代表近红外波段的反射率； ${\rho }_{Red}$ 代表红波段的反射率。NDVI能够很好地反映植被生长状况，取值大小与农作物生长茂盛程度相关，一般作物取值在0~1之间，约接近1作物生长越茂盛。

2.3.3. 遥感影像分类算法

本研究主要采用监督分类方法中的随机森林算法对种植结构进行识别。基本思路是在收集到的原始样本中，随机抽取一定量样本组成决策树，这些决策树具有样本和特征空间，抽取之后把样本放回原始样本中，使得每棵决策树内部都可能存在相同的样本和特征，实现每棵决策树都有对应的随机且独立数据集，以此增加算法的随机性，提高模型的泛化能力，之后将所有的决策树组合，对未知个体进行识别，最终对所有决策树分类的结果进行对比，选取个体的最优分类作为最终的输出结果。

一般抽取每个样本集的原则是不超过总样本数的三分之二，相对于每次抽中的样本而言，约有三分之一的样本未被抽中，利用这一部分未被抽中的样本估计训练的内部误差 [2] ：

$Q_i=\frac{M_i}{S}<\frac{2}{3}i=1,2,3$ (3)

式中： $i=1,2,3\cdots n$ 为抽取样本的次数， $Q_i$ 表示抽中样本占总样本的百分比， $M_i$ 表示每次抽取的样本数， 为原始样本数。然后根据样本建立各自对应的分类决策树，并将所有的决策树合并为一个森林，基于每棵决策树对未知个体所属类别进行判断，统计每一类作物在所有树的识别结果中所出现的频率，出现最多的作物类别就认定为该个体最终的分类结果 [2] ：

$p\left(x\right)=\underset{i=1}{\overset{K}{{\displaystyle \sum }}}{P}_i\left(x\right)$ (4)

式中：K为随机森林中树的个数，x为类别，P_i表示在决策树中该农作物类型在分类过程中所出现的比重。

2.3.4. 样本采集

本文采用实地勘测的方式采集样本点并使用目视解译法样本点进行扩充，对典型作物进行提取以及精度验证。采用亚米级手持GPS (ICE980RTK)进行研究区种植类型野外样本调查，仪器坐标定义为WGS-84坐标系，测量值采用地理坐标，测量时间为2020年7月3日~2020年7月17日，询问当地农民2020年上半年种植冬小麦的情况之后再采集样本点，每种类型各采样20个，共采集120个样本。根据收集到的120个野外样本的特征信息采用目视解译法在地图上进行扩充，共得到22,570个样本点，其中70%作为训练样本，30%作为验证样本。

2.3.5. 精度评价

使用验证样本对分类结果进行评价，通过混淆矩阵进行计算分析，获取提取结果的精度。混淆矩阵基于像元尺度进行计算，通过对比真实像元与影像对应像元位置的分类结果，利用相应指标将结果表示出来，主要的精度评价指标分为总体分类精度OA、Kappa系数、用户精度和制图精度。总体分类精度是指每一类验证样本所对应的分类结果的正确个数与所有实验验证样本个数之间的比值；Kappa系数是对一致性程度的评价判断指标，大于0.8表示分类结果较好。具体计算公式如下 [21] ：

$P_{0A}=\frac{{{\displaystyle \sum }}_{i=0}^n{x}_i}{N}$ (5)

$Kappa=\frac{P_{总体}-\frac{{{\displaystyle \sum }}_{i=0}^n{a}_i\times b_i}{N^2}}{1-\frac{{{\displaystyle \sum }}_{i=0}^n{a}_i\times b_i}{N^2}}$ (6)

式中： $P_{0A}$ 为总体分类精度； $x_i$ 为地物中被正确分类的个数； $a_i$ 为n类地物中的真实样本像元数； $b_i$ 为样本中被归为n类地物的像元数；n为分类结果总类别数；N为样本像元数。总体分类精度表示样本于分类结果的相似程度，Kappa系数的计算结果通常处于0~1之间。小于0.6则表示分类效果较差；位于0.6~0.8之间，分类精度较高；大于0.8表示分类结果与真实地物具有高度一致性。

3. 结果与讨论

本文通过随机森林模型分类法得到的6种典型作物种植结构分布如图2所示。

Figure 2. Spatial distribution map of six typical crops in Guangping County

图2. 广平县六种典型作物的空间分布图

从图中可以看出，该县冬小麦–夏玉米和夏玉米种植面积占比超过70%，其余农作物零星分布。冬小麦–夏玉米的种植分布范围最广，其中南阳堡镇、十里铺镇、南韩镇的种植比例较高。玉米和棉花在县域北部及南部分布较多。玉米绝大部分分布在平固店镇、东张孟镇；棉花主要分布在平固店镇南部、十里铺镇西部和胜营镇南部。花生、红薯和谷子的面积较为分散，在各镇均有少量种植。经统计，6种典型作物面积如表1所示。

Table 1. Area of various crops

表1. 各类作物的面积

如表所示，广平县内冬小麦–夏玉米、棉花、玉米、红薯、谷子和花生的种植面积分别为148.23 km²、9.34 km²、18.25 km²、2.31 km²、0.77 km²和3.98 km²。其中，种植比例分别为81.05%、5.11%、9.98%、1.26%、0.42%和2.18%，这与《2021年邯郸统计年鉴》上面的数据较为吻合，说明提取结果较为合理。分类作物混淆矩阵如表2所示。

由表可知，分类结果的总体精度为89.95%，Kappa系数为0.88，总体精度和Kappa系数均高于0.85，说明总体分类质量较精准。冬小麦–夏玉米、玉米和棉花制图精度都高于90%，冬小麦–夏玉米更是高达93.01%，说明这几种农作物的类别错分现象极低，分类结果基本正确。而制度精度稍低的红薯、谷子和花生由于种植面积较小验证样本总数少，导致的错分一个点对制图精度的干扰极大，但整体的用户精度都高于80%，说明其他类别之间相互错分的几率极低。

Table 2. Classification of crops classification confusion matrix

表2. 分类作物分类混淆矩阵

4. 结论

利用随机森林与哨兵2号遥感影像成功获取了2020年研究区内6种典型作物的空间分布图。与传统的方法相比，在提高精度的同时增加了作物的种类。冬小麦–夏玉米、棉花、玉米、红薯、谷子和花生提取之后的制图精度分别为94.95%、93.21%、91.76%、90.30%、88.92%、88.58%；用户精度分别为93.01%、91.5%、92.78%、89.85%、89.29%、87.70%，总体精度达到了89.95%，Kappa系数为0.88。各精度均满足要求，说明基于随机分类算法进行种植结构分类稳定性高，值得在生产实践中进行应用。经统计，冬小麦–夏玉米与玉米的面积最大，约占整个研究区总面积的70%左右，其余农作物呈零星分布。经统计，冬小麦–夏玉米、玉米、红薯、谷子、棉花和花生的种植面积分别为148.23 km²、9.34 km²、18.25 km²、2.31 km²、0.77 km²、3.98 km²和24.4503 km²。