陶万成, 张 颖, 谢茈萱, 王新盛, 董 镱,张明政 , 苏 伟*, 李佳雨, 轩 阜

1.中国农业大学土地科学与技术学院，北京 100083 2.农业农村部农业灾害遥感重点实验室，北京 100083

引 言

黑土是宝贵的自然资源，是肥力最高、最适宜农耕和最具生产潜力的土壤。黑土地保护性耕作模式是保护黑土地，促进农业可持续发展的重要措施。保护性耕作是指作物收获后地面上秸秆覆盖度大于30%的耕作和种植系统[1]，针对黑土地的保护我国东北地区正积极实施、推广保护性耕作模式[2]。秸秆还田不仅有利于资源合理利用，还能避免因秸秆露天焚烧造成的环境污染[3]，与传统耕作方式相比，长期实施秸秆还田有助于降低土壤容重，提高土壤孔隙度、土壤有机碳和土壤大团聚体百分比等[4]，对土壤保护与资源可持续利用起着重要作用。在亟需对黑土地实施保护措施的形势下，究竟哪些地块、有多少黑土地实施了秸秆覆盖，对黑土地保护相关科学研究和当地政府政策制定与实施具有重要意义。

传统作物秸秆覆盖区地面调查和识别的方法费时、费力、成本高，难以在区域范围内实施。遥感是快速获取区域范围内地表信息的技术手段，特别是欧空局(European Space Agency，ESA)发布的Sentinel-2卫星影像具有对秸秆纤维素和木质素敏感的短波红外波段(SWIR)，是快速、准确监测区域范围内黑土地上作物秸秆覆盖的有效数据源。已有许多研究者根据秸秆以2 100 nm处的强吸收谷与木质素、纤维素的高度相关性，利用多光谱影像、地面实测秸秆高光谱数据构建光谱指数展开秸秆覆盖度研究[5-7]。秸秆覆盖度研究一般仅利用作物收获后单一时相影像估算，无法表征不同秸秆覆盖类型的实际秸秆覆盖量，例如，高留茬秸秆覆盖从覆盖度的角度分析其覆盖度并不大，但实际上是4～50 cm的留茬直立在地表，其实际覆盖量并不小。因此，本研究进行秸秆覆盖区识别，是进一步开展地块尺度秸秆覆盖的基础，也是对秸秆覆盖类型精细分类(如覆盖还田、粉碎还田、高留茬等)的前提，通过遥感影像直接获取高精度的玉米秸秆覆盖区更具有实际意义。目前，利用遥感技术的秸秆覆盖区识别研究还较少，且没有综合考虑作物生长期和收获后的状态信息。本工作基于作物不同状态信息对玉米秸秆覆盖区开展精细、智能识别研究。

玉米是我国东北黑土地上主要种植的一种作物；吉林省四平市内有针对黑土地保护创建的“梨树模式”区域，广泛存在少耕、免耕等保护性耕作方式，又大量存在翻耕等传统耕作方式，非常适合开展秸秆覆盖区识别研究。由于玉米秸秆覆盖地块的反射率较高，且具有大面积连通的自然特性，在Sentinel-2A影像上表现为亮度偏高特征，呈现区域规则连续分布。根据秸秆覆盖大于30%时为保护性耕作的规定，遥感影像可视为秸秆覆盖区和背景两类。随机森林通过给定样本的学习训练形成分类规则，具有分析复杂地理信息系统分类特征的能力，适用于分类和高维数据的回归分析，同时具备更高的精度和稳健性。在作物倒伏识别[8]和土地覆盖分类[9]等方面，随机森林都具有很好的分类效果和处理速度，因此受到广泛关注。基于像素的随机森林算法分类结果中普遍会出现细小图斑。为克服不足，应用形态学滤波去除图斑。

鉴于此，本研究基于GEE云平台[10-11]，综合考虑玉米收获前、后的状态信息，结合研究区内2020年5月—11月份的时序Sentinel-2A影像创建的光谱指数和分位(quartile, QT)特征构建数据集，应用优化后的随机森林方法准确识别玉米秸秆覆盖区域；根据秸秆覆盖大范围连通的特性，利用连通域标定法优化分类结果，实现玉米秸秆覆盖区制图。

1 实验部分

1.1 研究区概况

研究区为吉林省四平市，地理范围为42°31′—44°09′N，123°17′—125°49′E，总面积约1.4万km2，地理位置如图1(a)所示。四平市处在黑土地和吉林省“黄金玉米带”核心地区，从2007年开始探索并实施以“秸秆覆盖、条带休耕”为主要内容的保护性耕作。开展区域范围内的玉米秸秆覆盖区准确识别，可应用于评估保护性耕作模式的实施程度，为监测保护性耕作推广，相关政策拟定提供可靠的数据支持。

图1 研究区位置与Sentinel-2A影像(日期：2020年11月4日)(a)与秸秆覆盖区照片(b)

1.2 数据源

1.2.1 地面调查数据

2020年11月3日—11日，在研究区进行了田间玉米秸秆覆盖实地调研。11月为研究区内玉米收获后秸秆覆盖的稳定时期，存在大量的样本地块。根据Sentinel-2A影像的最高空间分辨率(10 m)，在样本地块中随机选取10个分布均匀、大小为1 m×1 m的样方，采用拉绳法测量。秸秆调查时使用两根以0.1 m间隔作为标记的1 m刻度绳，垂直摆放在样方中，计算标记点与地面秸秆相交点数量与总标记数比值，统计10个样方的均值；依据保护性耕作的定义[1]对秸秆覆盖区进行识别，图1(b)中当统计的样方均值大于30%时为秸秆覆盖区，主要由机械收获后产生，包括低和高留茬；低于30%时为非覆盖区，主要是人工收获后产生的根茬。为了识别玉米秸秆覆盖区，采用实地调查和目视解译相结合的手段获取样本数据集，在剔除异常值后保留了646个样本：包括玉米秸秆覆盖、水稻秸秆覆盖、建筑用地、树林和水体样本，样本数量分别为200，103，118，137和88。

1.2.2 Sentinel-2A遥感影像

从Google Earth Engine(GEE)云平台获取时序Sentinel-2A影像，时间分辨率为5 d，光谱包括可见光、近红外、红边、水汽、卷云和短波红外波段，各波段分辨率分别为10，20和60 m。为避免低分辨率波段对分类任务的影响，舍弃空间分辨率为60 m的水汽、卷云和可见光中的气溶胶波段，选择B2(Blue)，B3(Green)，B4(Red)，B5(Red-edge1)，B6(Red-edge2)，B7(Red-edge3)，B8(Nir)，B8A(Red-edge4)，B11(SWIR1)和B12(SWIR2)较高空间分辨率波段。

利用GEE获取研究区内2020年的云占比小于10%的Sentinel-2A影像；同时选取20个实地调查的玉米秸秆覆盖样本点，统计年内归一化差值植被指数(normalized difference vegetation index, NDVI)和归一化差值秸秆指数(normalized difference residue index, NDRI)特征曲线和置信区间，如图2所示。S1为玉米出苗期，H为玉米收获期，S2为玉米秸秆覆盖期，横坐标为年积日(day of year, DOY)，图2(a)纵坐标为NDVI均值(NDVIM)和NDVIM+/-标准差(Std)，图2(b)纵坐标为NDRIM+/-Std。图2(a)中DOY为204时达到NDVIM最大值0.864。图2(b)中NDRI和相应时间区间的NDVI趋势相似，但值范围相差较大。为降低计算量，选取从5月13日到11月29日的182景云占比小于10%的Sentinel-2A作为影像源，充分合理地利用时间序列特征。根据研究需求，对获取影像进行波段筛选、拼接和裁剪等预处理。

图2 研究区2020年内玉米种植区时序NDVI和NDRI变化曲线

1.3 方法

基于时序Sentinel-2A的玉米秸秆覆盖区识别包括3个步骤。首先，基于Sentinel-2A时序影像计算光谱和指数特征，同时利用分位法获取分位QT特征，进而构建数据集；然后，对样本数据按7∶3随机划分为训练集和测试集，结合RF模型分类出玉米秸秆覆盖区域，而基于像素的RF分类结果中极易产生细小图斑，研究采用连通域标定法优化全局；最后，通过评价指标定量分析结果。具体技术流程如图3所示。

图3 玉米秸秆覆盖区识别技术流程图

1.4 数据集构建

1.4.1 光谱特征

对研究区内2020年11月4日获取的Sentinel-2A多光谱反射率影像进行统计分析得知，建筑物光谱特征较明显、容易识别，因此重点分析玉米秸秆覆盖区、水稻秸秆覆盖区、树林和水体样本光谱特征。统计样本内各地类波段反射率值并绘制箱线图，结果如图4(a—d)。

从图4(a,b)可以看出，玉米秸秆和水稻秸秆的光谱反射率区间范围为0.10～0.45，图4(c,d)水体和树林光谱反射率区间范围分别为0.00～1.05和0.00～0.30。通过对比发现，水体和树林相对秸秆反射率明显较低，易于分类。水稻秸秆和玉米秸秆的光谱反射率变化趋势和区间范围基本相同，B2—B8波段反射率呈上升趋势，B8A—B12波段呈下降趋势，仅利用单期遥感影像很难进行地类判别。鉴于此，考虑到不同作物的物候和结构特性，本研究使用7月22日Sentinel-2A影像提高水稻和玉米地类的分离度，即在作物生长季中玉米种植区的NDVI值大于水稻。为了判断该期影像是否有助于区分不同作物，利用样本统计绘制出玉米和水稻的光谱反射率箱线图，结果如图5(a,b)。图5(a)中玉米B2—B5波段反射率明显低于图5(b)水稻，B7—B8A反射率高于水稻，其他波段反射率相似，其原因可能与两种农作物的冠层结构相关，这些特征有利于分类。因此，光谱特征集是结合7月22日和11月4日的Sentinel-2A影像建立。

图4 不同地物类型光谱反射率特征(2020年11月4日Sentinel-2A影像)

图5 不同作物光谱反射率特征(2020年7月22日Sentinel-2A影像)

1.4.2 指数特征

为区分玉米秸秆覆盖区、水稻秸秆覆盖区、裸地和其他相似地类，分别构建了差值植被指数(NDVI)和秸杆指数(NDRI)。其中，NDVI可以较好地反映出不同农作物的生长状况，因此用于区别玉米和水稻种植区域。NDRI光谱指数与NDVI相似，同样是利用双波段进行计算，该光谱指数与玉米秸秆覆盖度显著相关[12]，其计算公式如式(1)

(1)

式(1)中，ρB12是短波红外2反射率，ρB4为红波段反射率。在玉米秸秆中含有大量的纤维素和木质素，反射率在2 100和2 300 nm附近有强吸收，即Sentinel-2A影像的B12波段反射率。使用Sentinel-2A的B12和B4波段反射率构建的NDRI可以最大限度地减少绿色植被的影响。综上分析，在光谱特征集的基础上，利用NDVI和NDRI建立指数特征集。

1.4.3 分位特征

分位QT特征的建立主要基于四分位数的思想。在统计学中，将所有数值按大小排列并分成四等份，处于三个分割点位置的数值为四分位数。同理，四分位数法可应用于时序遥感影像以构建QT特征，流程如图6所示。首先，通过GEE云平台获取时序影像集，图6(a)中为11月4日影像中蓝波段热度图，颜色越深反射率值越小，相反则反射率值越大；以中间红框内的像素为例，抽取该位置时间维度上的数据，如图6(b)；对获取的数据按照从小到大进行排序，V1—V5分别对应0%，25%，50%，75%和100%，依次选择对映分位上的反射率值，如图6(c)；最后，按照上述方法对影像中像素遍历处理获取光谱分位集，如图6(d)。

图6 分位QT特征构建示意图

QT特征还可以应用于时序指数特征集，与计算光谱分位集相似，首先获取时序指数特征集，然后应用QT特征提取指数分位集。QT特征集构建包括光谱分位集和指数分位集，这种提取特征的方式利用时序特性的同时避免了数据冗余，极大地降低了计算量。

1.5 随机森林秸秆覆盖区识别

随机森林监督分类算法是一种基于多棵CART决策树组合构成的集成学习方法。为了完全划分变量空间，利用自助抽取输入样本和节点随机分裂技术构建多棵决策树。通过单个决策树的多数投票策略获得预测结果，有效避免单个模型或参数组引起的误差，具有更高的精度和鲁棒性[13]。由于其分类的优越性，该分类器得到遥感界的广泛认可。

针对所设计的数据集，具有高维度特性，随机森林可以将每个维度的重要性输出到分类中，帮助特征选择以提高效率。因此，本研究使用嵌入在GEE中的随机森林算法对玉米秸秆覆盖区进行分类。根据分类结果的准确性调整决策树数量以期生成最佳的分类模型。

1.6 连通域标定和评价指标

研究区中秸秆覆盖地块一般具有大面积连通的自然特征，而本研究基于像素级别对玉米秸秆覆盖区进行分类，结果中会出现细碎的图斑，因此有必要借助连通域标定法对结果全局优化。连通域标定法示意图见图7，假定图7(a)为分类结果中某一区域，标签1表示为玉米秸秆覆盖区，标签0表示背景类，绿色区域为设定的4连通域滑动窗口。通过滑动窗口寻找标签为1的连通域，保留图7(b)左上角最大连通域(黄色区域)，删除右下角最小连通域(黄色区域)。如果图像中存在大量目标区域，则需标定所有的连通域，通过设置阈值对全局处理。

图7 连通域标定示意图

Kappa系数是衡量分类结果空间一致性的指标，明确揭示了分类结果的空间变化。总体精度(overall accuracy, OA)以预测正确与总体数量之间的比值反映分类结果的质量[14]。本研究选择上述指标来定量定性评价分类结果。

2 结果与讨论

基于时序Sentinel-2A影像设计的数据集，采用随机森林算法识别玉米秸秆覆盖区。为确定算法中最佳决策树数量，研究对比5，10，20，30和40决策树对分类模型的影响，具体见表1。从5到30整体呈上升趋势，30达到最大值并趋于稳定，当决策树数量设置为40时，评价精度略有下降。因此，实验统一标准，将决策树数量均设置为30。

表1 不同决策树的分类模型定量评价

2.1 基于不同特征集组合分类

在相同训练样区条件下，通过不同的特征集组合识别玉米秸秆覆盖区，定量评价结果见表2。以光谱特征集为基础，模型M5精度最高，Kappa/OA为97.41%/97.91%，比M1，M2和M4模型的Kappa/OA分别提高了5.84%/4.69%，4.52%/3.64%和1.29%/1.04%，表明相比于指数特征集，加入QT特征集可以明显改善模型精度；以QT特征集为基础，M4比M3模型提高1.95%/1.56%，表明与指数特征集相比，加入光谱特征集更有利于分类；以指数特征集为基础，M3比M2模型提高1.28%/1.04%，表明加入QT特征集优于光谱特征集。

表2 基于不同特征集的分类定量评价结果

本研究使用五种不同特征集组合模型，测试结果见图8，图8(a)和图9(a)分别为2020年11月4日和7月22日真彩色合成图，分类结果中蓝色和白色区域分别表示玉米秸秆覆盖区和背景。根据图8(b)—(f)对比发现，M5模型一定程度上降低细小图斑的产生，同时保留道路、农田等边缘细节信息，效果最好；M1和M2模型分类结果较差，易产生噪声，将道路等地类分为秸秆覆盖区；M3和M4模型分类结果相对较好，保留了部分边缘细节信息。实验结果与定量评价结果基本一致，QT特征的加入极大改善了分类可视化结果。

图8 基于不同特征集的分类可视化结果

2.2 基于不同时序尺度的特征分类

通过上述实验证明了QT特征的重要性，为了评估不同时序尺度生成的QT特征集对识别秸秆覆盖区的影响，本研究基于M5模型，对利用不同QT特征集的M5_1—M5_6模型定量分析，具体信息见表3。

根据表3，利用5月—11月影像创建的QT特征M5_6/M5模型分类效果最好。其中，M5_2比M5_1模型的kappa/OA低0.97%/0.77%，可能是10月为玉米收获期，秸秆地块中包括收获与未收获，收获的地块又包括秸秆成捆或散放在地里等，信息多样，识别困难；M5_3相比于M5_2模型精度低，分析认为9月份研究区受到台风影响，缺少质量好的影像。而M5_4，M5_5和M5_6精度相继提升表明长时序影像创建的QT特征集有利于识别秸秆覆盖。

表3 基于不同时间尺度特征分类定量评价结果

2.3 结合连通域标定法分类

在M5模型基础上，结合连通域标定法构建模型M6对全局优化，设置连通域阈值为30，M6分类结果的Kappa/OA为96.76%/97.36%。相比于M5模型精度有所降低，但仍高于其他模型，同样适用于玉米秸秆覆盖区识别。选取两组测试结果对比如图9和图10。

根据图9(c)，M5[图9(d)]和M6[图9(b)]相比于M5_1模型分类结果可以很好的保留细节信息，进一步证明QT特征设计的有效性；M6相比于M5模型在保留绝大部分细节信息的同时，抑制了细碎图斑的产生。图10(a)是夏季的子图块，其中浅绿色为水稻地块，图10(b)是11月初秸杆覆盖的子图块通过对比图10(c)M5_1，图10(d)M5和图10(e)M6分类结果发现，在仅利用11月影像的M5_1模型容易将水稻秸秆覆盖区分为玉米秸秆覆盖区，而利用5月—11月份影像的M5模型可很好地将两者区分开。此外，M5模型结果中有相对较少的细碎图斑，结合连通域标定的M6模型在优化细碎图斑后得到更理想的分类结果，见图10(e)。

图9 结合连通域标定分类结果1

图10 结合连通域标定分类结果2

通过实验对比M6模型可视化效果最好，定量评价结果仅次于M5模型，因此研究使用M6模型对整个研究区进行预测，得到的分类结果见图11。其中，左侧为整个研究区识别结果，右侧为研究区框选出的子区域识别结果，图11(a)主要包括秸秆和建筑，图11(b)主要包括秸秆，建筑，水域和树林。从不同区域识别结果可以看出秸秆覆盖区边缘信息保持完好，几乎没有细碎图斑，证明了所提出方法的有效性，适用于玉米秸秆覆盖区识别。

3 结 论

以吉林省四平市为例，基于GEE云平台，首先对时序Sentinel-2A影像预处理，利用光谱，指数和QT特征构建数据集。然后应用最佳决策树个数的随机森林分类器分析不同特征组合，不同时序尺度的QT特征和是否结合连通域标定对识别玉米秸秆覆盖区域的影响。最终采用可视化效果最好，精度较高的M6模型对研究区影像预测，主要研究结论如下：

(1)利用不同特征集组合的分类模型所设计的分类模型M5定量评价结果最优，通过模型对比表明加入QT特征能够有效提升精度，保留边缘细节信息，在一定程度上降低了噪声出现的概率，优于未加入QT特征的分类结果。

(2)基于不同时序尺度的QT和光谱特征集的分类模型，实验利用5月—11月时序影像内QT特征比较短时序内影像的QT特征分类效果好，表明长时序特征有利于玉米秸秆覆盖区识别，同时可以避免其他作物秸秆覆盖区域的干扰。

(3)在模型M5分类结果的基础上，结合连通域标定法的M6模型对全局连通域进行识别，通过设置阈值优化结果，实验结果表明该方法在保证较高玉米秸秆覆盖区域识别精度的同时去除了细碎图斑。

通过本研究提出M6模型可以准确的从Sentinel-2A影像中识别出玉米秸秆覆盖区，对当地监测保护性耕作实施和相关服务政策的拟定提供便利。后续研究主要有两点：① 将着重挖掘影像空间域纹理特征在玉米秸秆覆盖识别中的应用潜力。②尝试多源遥感影像融合进一步优化结果。