松嫩平原典型农区玉米秸秆覆盖度遥感估算

项小云,杜 嘉,赵博宇,周昊昊,宋开山

(中国科学院 东北地理与农业生态研究所，吉林 长春 130102)

0 引 言

农作物秸秆管理，是当今保护性耕作技术运用最为重要的一种保护措施[1]，秸秆覆盖有利于作物生产和环境的可持续性，可以有效减少土壤侵蚀、提高水分利用效率和土壤肥力[2]。秸秆覆盖度的估算是农作物秸秆管理的重要组成部分，秸秆覆盖度信息的获取，不仅可以了解保护性耕作的空间分布状况，还可以对于保护性耕作实施的进程及实施范围进行宏观的监测，提高政府监管的效率。因此，对农作物秸秆覆盖度进行估算能够为保护性耕作技术的实施提供数据支撑。

目前，秸秆覆盖度数据的获取主要依赖于人工现场数据收集和调查，使用这些方法进行大范围准确、系统和连续地获取数据是极其困难的。遥感技术可以系统的对秸秆覆盖度数据进行低成本有效获取，具有快速准确的特点，可以实现对秸秆覆盖度数据连续和大范围的获取[3]。Gausman等[4]较早开始利用 MSS 影像，根据光谱反射特征进行土壤与秸秆的区分。Daughtry[5]测量了玉米、大豆等农作物及土壤覆盖在400～2 400 nm波长范围内的地物光谱反射率，确认了使用纤维素吸收指数(Cellulose Absorption Index，CAI)进行秸秆覆盖识别的可行性。李志婷等[6]分析了Landsat-8 OLI各波段反射率及不同光谱指数与小麦秸秆覆盖度的关系，结果表明Landsat-8 OLI构建的光谱指数优于Landsat-5 TM构建的光谱指数。这些研究表明了遥感技术在秸秆覆盖的识别以及估算方面应用的可行性。

农作物秸秆的残茬和土壤有着相似的光谱，由于作物残留物具有与纤维素和木质素相关的性质，所以在2 100 nm附近形成独特的吸收特征[5]。这种独特的吸收特征成为基于光学遥感图像区分作物残余物与土壤的基础。因此，耕作指数的设计重点在于放大作物残留吸收特征的信号，同时抑制土壤及绿色植被等混淆的光谱信号[1]。以往研究中，诸多学者基于多光谱和高光谱遥感数据以及光谱辐射地面测量数据，构建不同的耕作指数方法模型对不同地区的作物秸秆覆盖度(Crop Residue Cover,CRC)进行估算。McNairn和Protz[7]使用Landsat-5 TM数据分析了CRC与归一化差异指数(Normalized Difference Index,NDI)之间的关系，结果表明：NDI与玉米秸秆覆盖度(MRC)相关性较好(R2=0.65～0.84)。Deventer等[8]分析了归一化差异耕作指数(Normalized Difference Tillage Index,NDTI)、简单耕作指数(Simple Tillage Index,STI)与CRC之间的关系，结果表明这些参数之间具有良好的关系。Zheng等[9-10]基于Landsat TM和ETM+数据，提取最小NDTI值(The minimum values of Normalized Difference Tillage Index,min NDTI)进行CRC估算，表明min NDTI和CRC之间关系较好(R2=0.89)。目前，基于遥感技术和耕作指数进行CRC的相关研究正趋于成熟，利用耕作指数进行CRC估算研究其精度的提高已经到达一定的高度，但是对于CRC估算精度的提高仍需进一步的摸索。

以往的CRC估算研究中，学者们多关注数据的光谱信息，忽视了对于空间信息的利用。纹理特征能够反映出图像亮度的空间变化情况，不同的地物类型之间纹理特征各不相同，能够较为直观地反映出影像的空间信息[11-12]。Li等[11]使用0.5 m分辨率的航空影像，进行了金银花的空间纹理特征提取并进行分类，结果表明了该方法分类的准确性。Wood等实地收集了193个采样点植被结构，使用两种不用数据影像提取纹理特征的信息进行植被结构的宏观测量，结果表明了使用纹理特征进行植被结构估算的可行性[13]。Najafi等人基于Landsat OLI数据，使用耕作指数、均值和纹理特征三种方法进行CRC的估算，精度和kappa系数分别为0.93和0.90、0.91和0.86、0.60和0.35，说明使用纹理特征方法进行提取和估算秸秆覆盖度具有可操作性[14]。

Zhou等[15]使用Quickbird数据进行叶面积指数(Leaf Area Index,LAI)的估算，评估了光谱指数、纹理特征参数和光谱指数与纹理特征组合方法三种不同的方法，结果表明：使用光谱指数和纹理特征组合方法可以显著提高LAI的估算精度(R2=0.84)。Berberoglu等[16]利用神经网络综合光谱和纹理信息进行地中海土地覆盖分类制图，表明组合方法可以提高分类精度。金秀良等[17]使用光谱指数和纹理特征组合方法进行秸秆覆盖度的相关研究，证明了运用组合进行秸秆覆盖研究的可行性。Ding等基于Sentinel-2数据构建光谱指数(Crop Residue Indices，CRIs)和纹理特征组合方法进行CRC的估算，三个MSI波段和三个纹理特征组合模型表现最佳(R2=0.69，RMSE=6.149)[18]。因此，光谱指数和纹理特征的组合方法在提高模型精度以及秸秆覆盖度估算方面具有应用价值。

本文以Landsat-8 OLI遥感影像为数据源，综合利用耕作指数、纹理特征和组合方法构建基于偏最小二乘回归方法的秸秆覆盖度模型，对松嫩平原部分地区的玉米秸秆覆盖度进行估算，并分析研究区秸秆覆盖度的空间分布特征，研究结果将为农业管理部门推动保护性耕作提供真实可靠的数据支撑。

1 研究区与研究方法

1.1 研究区概况

研究区位于东北三大平原之一的松嫩平原，是国家重要的粮食储备基地，地跨吉林与黑龙江两省，是我国主要农作物种植区之一，玉米种植面积占耕地面积的80%[19-20]。研究区气候属于温带大陆性季风气候，降水主要集中于夏季，雨热同期有利于农作物的生长。研究区选择地理位置在松嫩平原东部地区，行政区划为吉林省长春市至黑龙江省哈尔滨市之间，为Landsat-8三景影像(118/28、118/29和118/30)覆盖的区域(图1)。

图1 研究区(Landsat-8影像，标准假彩色band5/band4/band3)Fig.1 Study area(false color composite Landsat-8 image R/G/B.band5/band4/band3)

1.2 数据获取

所选用的Landsat 影像来源于美国地质调查局USGS 网站(https://espa.cr.usgs.gov/index/)，通过查询地理位置匹配关系，筛选影像上基本无云或者云量较少的有效3景影像，时间为2019年2月3日。本研究选取的3景影像部分区域有少量云层覆盖，将该部分区域进行去云处理，以消除云层对反演产生的影响。并且3景影像时间为同一天，不涉及不同图像间的光谱匹配问题，因此直接对3景影像进行拼接合并。然后，我们利用2018年玉米种植空间分布数据对遥感影像进行裁剪，以达到和地面实测匹配的要求。

2019年3月27日至4月4日，对影像覆盖区玉米秸秆覆盖度进行野外调查，使用样线法[21]对秸秆覆盖度进行粗略估算。在充分考虑道路可达性的前提下，采样点选择距离道路和村庄较远的农田区域，并且选取采样的农田四周没有树木和电线杆等其他地物，避免其他地物产生的干扰。实地使用了50 m长的绳子，将其间隔1 m分成50个部分，利用红绳标记，开始进行测量。在每个采样区域，用两根绳子沿行对角线进行垂直拉伸，并计算与作物秸秆相交的标记数，用以计算该采样点的秸秆覆盖度。此外，使用GPS记录每个采样点的位置，获取采样点照片，完成信息的采集。调查共采集玉米采样点34个，将数据随机分为24个校准数据集和10个验证数据集两部分。

针对影像和实测数据间隔时间对研究产生的影响，主要有以下考虑。近几年来，东北三省地区严格贯彻落实国家环境保护部下发的秸秆禁烧政策，秸秆焚烧状况得到较为有力的控制。因此，在影像和野外实地测量的时间间隔中秸秆覆盖度变化的几率很小。

1.3 研究方法

利用耕作指数(Tillage Index,TI)、纹理特征和组合方法构建基于偏最小二乘回归方法(PLSR)的玉米秸秆覆盖度模型，对松嫩平原的部分地区进行玉米秸秆覆盖度MRC估算。同时，利用决定系数(The coefficient of determination，R2)和均方根误差(Root Mean Square Error，RMSE)进行统计与验证的分析。

1.3.1 耕作指数。耕作指数就是指基于玉米秸秆独特的光谱吸收特征，从而构建的用来进行秸秆覆盖度估算的一种光谱指数，是目前研究秸秆覆盖度估算的常用方法之一。在遥感中并没有针对该类光谱指数进行准确的定义，有少量相关的中文文献仅是使用光谱指数这一广泛的定义进行描述，本文使用的耕作指数是基于众多相关研究中的英文文献对此的描述进行直译[1,6-17]。

本研究中使用了7种耕作指数，分别为归一化差异耕作指数(NDTI)、简单耕作指数(STI)、归一化差异指数(NDI7、NDI5)、短波红色外归一化差异指数(SRNDI)、归一化差异衰老植被指数(NDSVI)、改良耕作指数(MCRC)，这7种指数是目前秸秆覆盖研究中使用较为成熟的几种(表1)。

表1 耕作指数及计算公式简介Table 1 Introduction of tillage indices and calculation formula

1.3.2 纹理特征。纹理特征指通过一定的图像处理技术提取出纹理特征参数，从而获得纹理的定量或定性描述的处理过程[24]。纹理分析在遥感图像处理中有着重要的作用，灰度共生矩阵(Gray Level Co-occurrence Matrix,GLCM)是纹理统计分析中最常用的方法之一[25-26]，使用灰度共生矩阵的目的就是用其辅助遥感影像纹理分类。本研究中选择8个纹理特征指标，包括均值(mean)、方差(variance)、同质性(homogeneity)、对比度(contrast)、相异度(dissimilarity)、熵(entropy)、二阶矩(second moment)和相关性(correlation)。这些纹理特征在进行地表生物量估算、植被结构测量以及土地覆盖分类中使用较为成熟，均证实了其在地表估算和识别中应用的可行性[13,16,27]。

通过ENVI5.3软件平台对Band2-Band7分别进行8种纹理特征值的提取。Band2-Band7波长范围在0.45～2.3 μm之间包括可见光、近红外和短波红外波段可用于秸秆覆盖的研究，其余波段不涉及该方面的应用。窗口的大小选择是进行纹理特征提取的关键，窗口的大小与值的范围有着直接的关系。使用较小的窗口会在一定程度上扩大窗口内的差异，但是可以保留较高的空间分辨率，而使用较大的窗口在提取过程中会使纹理变化过度平滑而导致提取效率偏低[28]。在考虑提取效率和分辨率的基础上，3×3像元为纹理特征提取的最小窗口，因此本研究选择窗口大小为3×3像元进行计算提取。这些统计指标可以在反映地物空间特征差异的同时从不同的角度来量化影像局部的纹理结构。

1.3.3 偏最小二乘法。偏最小二乘法是从实际应用出发的一种多元统计数据分析的新方法，这种方法主要用于多因变量与多自变量之间的线性回归模型构建[29]。这种方法的优势在于，当观察数量小于变量数量时会出现多重共线性的问题，对模型的构建产生较大影响，偏最小二乘法的使用可以消除由此产生的不利影响。本研究涉及耕作指数和纹理特征指标过多，并且各波段组合计算之间存在一定相关性问题。通过前期进行主成分分析以及计算其自变量间的关系，发现各自变量间存在不同程度的共线性问题。因此，采用偏最小二乘法就是为了消除多种自变量之间的共线性问题。

PLSR模型是多个线性回归模型的扩展。当需要从大量自变量中预测一组因变量时，PLSR模型具有特别的优势，同时可以有效地解决多自变量中的强共线性和噪声问题。在最简单的形式中，线性模型指定了从属(响应)变量y和一组预测变量(x)变量之间的(线性)关系[30]：

y=b0+b1x1+b2x2+b3x3+…+bpxp

式中：b0是回归系数的截距；bi值是数据计算从变量1到p的回归系数。

1.3.4 统计与验证分析。相关性是分析两个或多个变量之间亲密程度的重要指标，Pearson相关系数通常用于衡量两个变量之间是否存在线性关系，决定系数是相关系数的平方。本研究中使用决定系数(R2)来分析MRC与耕作指数和纹理特征之间的关系。研究通过比较R2值的差异来确定模型的性能。一般来说，较高的R2值代表预测MRC模型的精度和准确度更高。

均方根误差表示预测值与测量值之间的偏差程度，可以很好地反映回归模型估计MRC的准确性。一般来说，RMSE值越小，测量值与预测值之间的误差越小，回归模型估计的准确度越高[31]。

2 结果与讨论

2.1 玉米秸秆覆盖度(MRC)与耕作指数(TI)的关系

不同耕作指数与MRC之间分别存在着一定的回归关系(图2)。其中NDTI、STI与MRC相关性最好，R2均为0.86，RMSE分别为3.69、3.78(表2)。NDI7、NDI5与SRNDI与MRC相关性在0.5～0.3之间，NDSVI、MCRC与MRC之间相关性较不明显。TI与MRC相关性由高到低依次为NDTI、STI、NDI7、NDI5、SRNDI、NDSVI、MCRC。这些结果表明，NDTI、STI、NDI7、NDI5、SRNDI可以用于MRC的估算研究。

表2 玉米秸秆覆盖度(MRC)与耕作指数(TI)的关系Table 2 Relationship between maize residue cover(MRC) and tillage index (TI)

图2 归一化差异耕作指数(NDTI)、简单耕作指数(STI)与玉米秸杆覆盖度(MRC)的关系Fig.2 Relationship between normalized difference tillage index(NDTI),simple tillage index(STI) and maize residue cover(MRC)

研究所选取的7种耕作指数，除了NDSVI和MCRC以外，其余指数均与MRC具有较好的关系。Daughtry等[5]研究表明，在干作物残留物的反射光谱中2 100 nm附近的吸收特征十分明显，这是在土壤光谱中不存在的，这可能与作物残留物中的木质素和纤维素有关。因此，2 100 nm附近的光谱信息对于识别农作物秸秆至关重要。NDTI、STI、NDI7、NDI5、SRNDI等5种耕作指数均涉及Band7(2 100～2 300 nm)，所以这4种指数与MRC相关性较好。

Daughtry等[32]研究表明作物残留物中水分含量会降低所有波段的反射率，在1 450 nm和1 960 nm附近形成两个吸收带(波长>1 300 nm)，最高光谱是烘干残留物，最低光谱为水饱和残留物。所以，对于波长>1 300 nm，位于1 450 nm和1 960 nm附近的反射光谱信息提取用于含水的作物残留物识别具有重要意义。NDTI和STI的计算是使用Landsat-8 OLI数据的Band6(1 560～1 660 nm)和Band7(2 100～2 300 nm)，因此，STI和NDTI与MRC相关性较好(表3)，R2均为0.86。Wanjura等[33]使用多波段辐射计测量反射率，结果表明波段MMR 4(760～900 nm)在风化期间表现出最大的相对反射率变化，换言之，波段MMR 4与风化期作物残留物具有显著相关关系。本研究进行MRC估算所使用影像时期为2月份，农作物在上一年11月末前全部收割完毕，此时秸秆正处于作物残留物风化期。NDI7值是由Landsat-8 OLI数据的Band5(845～885 nm)和Band7(2 100～2 300 nm)计算得出，因此NDI7与MRC之间具有良好的相关性(R2=0.437)。Daughtry等[32]研究认为，基于2 100 nm附近的纤维素和木质素吸收特征可以用来很好的区分作物秸秆。这或许可以解释在本研究中SRNDI与MRC的相关性要高于NDSVI。因此，利用NDTI、STI、NDI7、SRNDI用于MRC的估算研究较为合理。

2.2 玉米秸秆覆盖度(MRC)与纹理特征的关系

研究选取了8个纹理特征指标(均值、方差、同质性、对比度、相异度、熵、二阶矩和相关性)，分别计算了Band2-Band7共6个波段的八种纹理特征的48个指标与MRC之间的相关关系。结果显示8个纹理特征指标与MRC之间的相关性较好。B5mean与MRC之间的相关性最高(R2=0.384)，与B4variance的相关性最低(R2=0.212)。8个纹理特征指标的R2值范围从0.384到0.212，且纹理特征与MRC回归关系由高到低依次为B5mean、B2mean、B4entropy、B4second moment、B4contrast、B4dissimilarity、B5correlation、B4variance(表3)。使用野外实测的覆盖数据验证MRC回归方程估算的准确性，比较了预测MRC值和测量值，RMSE结果在8.03%～17.59%之间。因此，根据统计结果综合考虑B5mean、B2mean、B4entropy可以用于MRC模型的估算。

表3 玉米秸秆覆盖度(MRC)与纹理特征的关系Table 3 Relationship between maize residue cover(MRC)and texture feature indicators

本研究中，8种纹理特征与MRC相关性相对较高。其中，B5mean、B2mean、B4entropy三种纹理特征与MRC相关性相对较高。剩余五种纹理特征与MRC之间的相关关系并不理想，这可能是由于农作物不同的秸秆覆盖方式带来的影响。通过实地调查和测量发现玉米的秸秆覆盖方式并不相同，有的区域在农作物收割完成后，将秸秆进行充分粉碎后均匀的铺在地表，地表纹理特征较为平整；有的区域在大型农作物收割机收获后，直接将秸秆残留在作物种植间隙之间，地表纹理特征呈现条带状；少数区域由于人工收割，秸秆随意的散落在地表覆盖，地表纹理特征较为凌乱。纹理特征正是对秸秆空间特征(空间信息)的一种表述[34]，由于秸秆覆盖方式的多样，使得纹理特征较为复杂，其中干扰的因素较多，从而没有得到较为理想的精度。

2.3 利用偏最小二乘回归(PLSR)进行MRC估算

研究使用TI、纹理特征以及组合方法三种方法，利用PLSR对MRC进行估算。首先，分别使用TI和纹理特征单独进行MRC模型的构建。对于TI方法来说，根据TI相关关系的高低进行选择。由于NDI5与NDI7均为归一化差异指数，且NDI7与MRC的相关性要高于NDI5，因此选择NDI7作为归一化差异指数指标。最终，选取NDTI、STI、NDI7、SRNDI四种指数进行MRC估算模型的构建。如表4结果所示，基于四种TIs的模型R2值为0.857，基于所有TIs的模型R2值为0.86。因此，向模型添加NDI5、NDSVI和MCRC并没有明显提高MRC估算准确性。对于纹理特征方法，根据纹理特征与MRC相关关系由高到低依次选择B5mean、B2mean、B4entropy三个指标利用PLSR进行MRC估算，模型R2值为0.185。

基于8种纹理特征指标(B5mean、B2mean、B4entropy、B4second moment、B4contrast、B4dissimilarity、B5corrleation、B4variance)的模型R2值为0.239，结果表明进一步向模型中添加指标可以提高模型的准确性，但提高效果并不显著。Li等[35]研究表明，使用PLSR进行建模时，只有少数几个特征便可以满足提取与区分基本信息，添加大量指标并不会显著提高模型估算的准确性，这与本研究结果相一致。当模型中只涉及4种耕作指数(NDTI、STI、NDI7、SRNDI)时，MRC模型估算的精度已经满足研究需求，进而添加耕作指数对于模型估算的准确性并没有明显提高，R2仅提高0.003。但是添加3种纹理特征进入模型，对模型的精度有所提高，这也进一步表明了耕作指数和纹理特征的组合是可以提高模型估算的准确性。

研究进一步尝试利用TI和纹理特征组合方法进行MRC估算模型构建。首先，选取4种TIs(NDTI、STI、NDI7、SRNDI)和3种纹理特征(B5mean、B2mean和B4entropy)指标进行组合，利用PLSR进行MRC估算模型构建，组合模型的R2值为0.849(表4)。进而选取所有的TIs和8种纹理特征指标(B5mean、B2mean、B4entropy、B4second moment、B4contrast、B4dissimilarity、B5corrleation、B4variance)构建模型，组合模型的R2值为0.847。研究发现，TIs和纹理特征的组合对于MRC估算模型的准确性并没有提高反而呈现降低的趋势。但是根据先前学者研究表明，TIs和纹理特征的组合可以显著提高模型精度[14-15]。因此，研究对于纹理特征指标的添加和组合进行了重新的选择，选择的指标依据金秀良等人先前在秸秆纹理特征的研究中梳理出的相关性较高的纹理特征。进而选取4种TIs(NDTI、STI、NDI7、SRNDI)和3种纹理特征(B3mean、B4mean和B5mean)指标进行组合，模型的R2值为0.881；所有的TIs和8种纹理特征指标(B2mean、B3mean、B3homogeneity、B3dissimilarity、B3entropy、B3second moment、B3correlation、B6mean)构建模型，模型的R2值为0.907。此时，TIs和纹理特征的组合可以提高模型精度。所以，使用组合方法的MRC估算精度高于单独的TIs和纹理特征方法，组合方法可用于MRC估算研究。

表4 基于偏最小二乘回归估算MRC三种方法的比较Table 4 Comparison of three methods for estimating MRC based on partial least squares regression

研究对耕作指数模型、纹理特征模型和基于两种组合方法构建的模型进行验证，我们使用野外实测数据和模型反演的MRC进行对比分析。如图3所示，除纹理特征验证R2较低(R2=0.510 3)，其余三种模型R2均大于0.7以上，可以满足研究的需求。

2.4 玉米秸秆覆盖度(MRC)空间分布

利用基于所有TIs和8个纹理特征指标的组合方法估算MRC的空间分布，整个研究区的MRC范围在8%～95%之间(图4)。

图4 松嫩平原MRC的空间分布Fig.4 Spatial distribution of MRC in the Songnen Plain

研究区内MRC平均覆盖度为45%，从行政区划来看，MRC较高地区集中在长春市和绥化市。长春市的九台市、榆树市及德惠市MRC较高，平均覆盖度分别为60.2%、52.4%、48.6%；绥化市的兰西县和肇东市MRC普遍较高，平均覆盖度分别为60.4%、58.1%。其余区域为哈尔滨宾县及四平市梨树县的MRC较周边区域覆盖度高。这主要是由于国家相关政策的实施，这些区域多为国家农业部建设的保护性耕作技术的示范基地以及保护性耕作技术实施较为广泛的区域[36]。部分地区由于保护性耕作部署方式不同，会在农作物收获之后进行秸秆深翻还田操作，在哈尔滨以北的区域在农业生产中存在着收割后随即进行秸秆深翻的农业措施，这也可能是导致部分地区秸秆覆盖度估算结果较低的原因[37,38]。

3 结论

NDTI、STI、NDI7、SRNDI与MRC相关性较好，其中NDTI、STI与MRC相关性最高，R2均为0.86。纹理特征B5mean与MRC之间的相关性最高(R2=0.384)。基于组合方法的MRC估算模型精度最高，7种耕作指数+8种纹理特征的组合模型精度R2可以达到0.907。

不同的纹理特征组合方法对于模型估算的精度影响很大。第一组纹理特征(B5mean、B2mean和B4entropy)指标与4种指数进行组合时，组合模型的R2值为0.849，纹理特征的加入反而略微降低了模型的精度。第二组纹理特征(B3mean，B4mean和B5mean)指标进行组合，组合模型的R2值为0.881，3种纹理特征的加入提高了模型估算的精度。7种耕作指数+8种纹理特征的组合方法中，第一组纹理特征(B5mean、B2mean、B4entropy、B4second moment、B4contrast、B4dissimilarity、B5corrleation、B4variance)加入，模型的R2值为0.847。第二组纹理特征(B2mean,B3mean,B3homogeneity,B3dissimilarity,B3entropy,B3second moment,B3correlation,B6mean)组合构建模型，模型的R2值为0.907。第二组选择的纹理特征相较第一组而言为MRC估算提供了更多的有用信息。

基于TIs和纹理特征的MRC估算方法的准确较高，这为区域尺度的MRC估算提供了一种便捷有效的新方法。鉴于本研究仅是基于三景影像的光谱信息和纹理特征对该区域内玉米秸秆覆盖度进行估算，今后的研究中应该获取更多Landsat-8 OLI影像的光谱信息和纹理特征，并且对模型指标的选择进行进一步研究，用以验证本研究结果对不同作物和环境的适用性。