基于无人机多光谱影像的完熟期玉米倒伏面积提取

张新乐，官海翔，刘焕军,2，孟祥添，杨昊轩，叶 强，于 微，张汉松※

(1.东北农业大学公共管理与法学院，哈尔滨150030：2.中国科学院东北地理与农业生态研究所，长春130012)

0 引 言

中国是一个农业大国，也是世界上受风灾影响较为严重的国家[1]。风灾倒伏面积的准确提取，能够为灾后的农业生产、政府决策以及农业保险理赔，提供数据和技术支持[2-3]。传统的作物风灾倒伏信息是通过田间实测的方式获取，工作量比较大；卫星遥感技术虽然能够进行全方面、大范围地提取倒伏面积信息，但时效性相对较差，且影像易受云层干扰[4-5]，而且一般卫星影像的空间分辨率难以满足地块尺度下的倒伏面积信息精准提取。近些年来，无人机遥感凭借精度高、受地形约束小、成本低、操作简便等优势，迅速占领了遥感市场，并逐渐成为农业灾害信息获取的重要方式[6-7]。与卫星遥感相比，无人机遥感不仅具有更高的时间和空间分辨率，而且能够获取田块尺度作物的空间变异信息，能反映地块内部不同位置和生长期作物的长势差异[8-11]。因此，无人机遥感技术非常适用于田块尺度的作物倒伏面积提取。

目前国内外研究主要使用无人机图像或卫星遥感影像，选取对倒伏作物识别度较高的特征对作物倒伏面积信息进行提取。Liu 等利用无人机获取的彩色图像和热红外数据，构建了倒伏水稻的识别模型，并对模型的识别精度进行了验证[12]。李广等对无人机彩色图谱进行了二次低通滤波分析，并将两种单特征的线性拟合作为特征集，由此提取了小麦倒伏的空间分布[13]。李宗南等利用Worldview-2 多光谱影像，分析了灌浆期倒伏玉米地块的纹理和光谱特征，基于优选的3 波段纹理和光谱特征估算了倒伏玉米面积[14]。Yang 等将无人机图像、纹理特征和数字表面模型进行了不同组合，并计算了不同组合特征中各类作物形态的分类精度，选取了彩色图像和数字表面模型组合特征作为倒伏水稻识别的指标[15]。李宗南等将研究区内的作物形态划分为倒伏玉米与未倒伏玉米两类，通过对该两类作物形态的纹理特征对比，选取了均值纹理特征作为倒伏玉米识别的指标[16]。Han 等使用无人机获取了倒伏玉米的多光谱和可见光图像，并提取纹理、冠层结构、植被指数等参数，以其为输入量构建了2 种倒伏提取的Logistic 模型，并验证了模型的精度[17]。Chu 等利用超空间分辨率的无人机影像构建了玉米的冠层高度模型，通过设立高度阈值，实现了倒伏与未倒伏玉米的分类[18]。毛智慧等通过无人机倾斜测量的方式获取了倒伏作物的数字表面模型（DSM），并结合无人机图像的色彩特征，进行了玉米倒伏信息提取[19]。刘良云等基于倒伏前后的两景TM 影像，利用NDVI 值监测了小麦抽穗初期的倒伏发生程度[20]。王立志等利用多时相的多光谱卫星影像，分析了抽雄期玉米倒伏前后的多种植被指数变化，采用RVI差值监测了倒伏受灾范围[21]。但以上研究主要集中于作物生长旺盛期的倒伏信息提取，而未对完熟期的作物倒伏提取进行研究，并且所使用的倒伏提取特征主要以单类纹理、光谱与植被指数特征或不同类型特征间的组合为主，弱化了纹理特征对于倒伏作物的空间域形态信息的反映，也未对纹理特征提取倒伏作物的最优维数进行探讨，这在一定程度上造成了倒伏识别特征的信息冗余和分类结果的过拟合效应。

正常生长的玉米在完熟期叶绿素含量较低，其叶片呈现淡黄色，与倒伏玉米叶片的颜色很相似，两者的反射光谱特征差异较小[22]。而处于生长发育旺盛期的玉米，其倒伏前后的光谱特征差异较大，因此基于生长旺盛期的玉米倒伏识别特征并不适用于完熟期玉米倒伏面积识别。本文利用5 波段的无人机多光谱数据，分析研究区内不同作物形态的光谱、植被指数及纹理特征，优选分类误差最低的特征组合，提取完熟期玉米地块内的倒伏面积，定量分析和评价倒伏面积提取精度，以期寻求最佳的完熟期玉米倒伏面积识别的特征组合和方法，为农业灾损的精准评估提供技术参考。

1 材料与方法

1.1 研究区概况

友谊农场位于黑龙江省双鸭山市友谊县西部，隶属于黑龙江农垦总局红兴隆管理局，研究区位于该农场西北部，常年种植玉米，占地面积为172 975 m2，地理坐标为131.6°E，46.7°N，具体位置如图1所示。此地块玉米正处于完熟期，并存在4 类典型作物形态，主要包括：长势较差、叶片呈现绿色的未倒伏A 类玉米；长势较好、叶片呈淡黄色的未倒伏B 类玉米；叶片呈现淡黄色的倒伏玉米；辐射亮度较低、分布零散的阴影区域。该地区位于三江平原地区，气候属于温带大陆性季风气候，雨热同期，冬季干燥，夏季温热多雨。平均年降雨量在540 mm 左右，年日照时数约为2 500 h，有效积温达2 300 ℃，海拔高度在64～70 m 之间，地势平坦，土壤类型为暗棕壤，肥力较高。主要种植玉米和水稻，一年一熟。

1.2 数据获取

本实验采用大疆600 M Pro 六旋翼式飞行器作为遥感平台，该平台质量为10 kg，最大可承受风速8 m/s，能够在-10～40 ℃的温度环境下正常工作，最大承重上限5.5 kg，最高续航时间是38 min。根据飞行器的荷载重量，搭载MicaSense RedEdge™3 多光谱相机，该相机有红、绿、蓝、红边、近红外5 个波段，配备光强传感器和一块白板。在飞行过程中，光强传感器可以校正太阳光线的变化对影像所产生的影响，白板则具有固定的反射率信息，利用白板对影像进行校正即可生成反射率图像。数据采集时间为2018 年9 月12 日，天气晴朗无云，风速小于1 级，适合飞行。为保证影像的完整和准确，通过地面站设置飞行的航向重叠为80%，旁向重叠为79%，主航线角度99°，边距-1.8 m，航线9 条，飞行高度为110 m，实际高度109.5 m，图像的空间分辨率为0.074 m，航拍完成后采用自动返航的方式降落。

图1 研究区位置及验证点分布Fig.1Locationandverificationpointdistributionofresearcharea

1.3 数据处理

首先对无人机采集的照片进行筛选，去除姿态角度异常、成像存在问题的图像。将筛选后的图像和POS 数据输入Pix4D mapper 软件中，根据相机的配置调整处理参数。运行软件后，应用会自动生成连接点并与POS 数据参与空三运算，由此得到每一张航拍影像的准确外方位元素和加密点的坐标。经点云加密后，软件会自动生成数字表面模型，以此作为数字微分纠正的依据对影像进行正射校正[23-24]。

1.4 研究方法

1.4.1 研究技术路线

研究首先基于无人机多光谱特征图谱分别选取4类作物形态的统计区域，提取并分析每一类统计区域的光谱、植被指数与纹理特征。然后经过特征筛选后分别构建适于倒伏提取的光谱、植被指数、单类纹理与多类纹理特征，利用以上4种特征组合分别结合最大似然法进行分类得到4种倒伏玉米面积的提取结果。最后验证每种结果的Kappa系数和面积误差，分析倒伏玉米面积提取精度，优选倒伏玉米面积提取的最佳方法。技术路线如图2所示。

1.4.2 特征统计样区选择

根据目视解译，依照玉米叶片的不同颜色将未倒伏玉米分为未倒伏A类和未倒伏B类；垂直特征不明显，垄向性较差的玉米归为倒伏玉米；黑色区域归为阴影。从多光谱影像上均匀选取未倒伏A类2 113.76 m2、未倒伏B类2 189 m2、倒伏区域2 164.07 m2，3类作物形态分别选5处，由于图上阴影区域分布零散且面积较小，因此阴影区域选择50处，共计253 m2。4类统计样区的选择应均质且不能包含其他类别。其中，未倒伏A类主要为地块东部受大风影响严重，叶表呈翠绿色，长势较差的未倒伏区域；未倒伏B 类为地块西部受大风影响较小，叶片表现为淡黄色特征，长势较好的未倒伏区域；倒伏区域为地块东部边缘，受大风影响而发生倒伏的区域；阴影区多为未倒伏玉米在倒伏玉米区域上的投影以及南侧防护林的阴影。

图2 倒伏玉米面积提取技术路线图Fig.2 Technology roadmap of maize lodging area

1.4.3 光谱特征选择

分别计算4 类特征统计样区中各个波段的反射率均值，对应各波段波长的中心位置，绘制每个类别的光谱反射率曲线，统计4 类区域在各波段上反射率的平均差异。虽然本研究使用的多光谱数据具有定量的光谱信息，但不同类别作物形态之间的反射率差异过小仍会产生较多的错分像元[25-27]，进而影响倒伏面积提取的精度，因此为增强各类区域之间的反射率差异，减少错分像元数量，本研究根据计算结果，选取平均差异相对较大的蓝、绿、红3个波段作为倒伏面积提取的光谱特征指标。

1.4.4 植被指数法

归一化植被指数（normalized difference vegetation index,NDVI）通过非线性拉伸的方式增强了红波段与近红外波段反射率之间的差异，是反映农作物生长状况和植被覆盖密度的重要指标[28-29]；而红边归一化植被指数（red edge normalized difference vegetation index, NDVIR-edge）则更能够反映叶冠层的细小变化、植被片断和衰老等信息，这些信息通过红波段与红边波段反射率的非线性拉伸而得以体现[30-31]。NDVI与NDVIR-edge的计算公式：

式中BRed为无人机多光谱影像中的红波段；BNir为近红外波段；B717为红边波段，其中心波长为717 nm。

基于无人机获取的多光谱影像分别提取NDVI 和NDVIR-edge，提取每类特征统计样区中的植被指数，分别统计2种植被指数中每类作物形态的均值。

1.4.5 纹理特征选择

对影像进行灰度共生矩阵的纹理滤波，得到各波段的均值、方差、协同性、对比度、相异性、信息熵、二阶矩、相关性等8类，共40项纹理特性，并统计各类区域中每项纹理特征的均值[32-33]。以统计结果作为分析的原始数据，分别进行基于单类和多类纹理特征的倒伏面积提取特征选择。

1）单类纹理特征选择：在只选取单类纹理特征作为倒伏面积提取指标的条件下，为使提取的倒伏面积更加准确，实验将基于8 类纹理特征，分别建立4 类作物形态的训练样区，利用最大似然法进行分类，根据分类结果选取适于区分4 类作物形态，倒伏面积提取精确的纹理特征。通过反复试验，最终选取均值纹理特征。

2）多类纹理特征选择：为了尽可能地降低纹理特征的维数，去除维数冗余所产生噪声信息，本文首先选取均值、对比度2类纹理描述，作为倒伏玉米面积提取的基础特征池，利用该特征池结合最大似然法进行分类，均匀选取倒伏和未倒伏玉米验证样本各50 个，计算提取结果的Kappa系数。然后，以该特征池为基础，以一类纹理向量为步长，对该特征集进行维数叠加，验证每次纹理向量叠加后分类结果的Kappa 系数，从而确定最优的分类维度，验证过程如图3 所示。根据图3 的结果，本研究最终确定多类纹理特征组合的最优维数为15 维，共计3 类纹理特征。

图3 不同维数下的Kappa系数变化Fig.3 Variation of Kappa coefficients at different dimensions

在确定用于分类的纹理特征之前，本研究将先对划分的未倒伏A类、未倒伏B类、倒伏以及阴影区的各项纹理特征进行对比分析，并统计4 类作物形态之间各类纹理特征的相对差异度。为了使所选纹理特征保有最大类间差异和最低类内离散度，进而提高分类特征的鲁棒性，本文选取了每2 类区域对比之下差异度较大的4 类纹理特征，分别为均值、对比度、二阶矩和协同性纹理特征。任选3 类特征进行组合，验证每种特征组合下倒伏玉米提取结果的Kappa 系数，选出Kappa 系数最高的纹理特征组合。经反复试验，最后确定用于分类的纹理特征分别为均值、对比度和协同性。

1.4.6 倒伏面积提取和精度验证

分别对NDVI、NDVIR-edge、单类纹理特征、多类纹理特征、光谱特征，共5 类特征集进行分类，由此获得每一类特征集下的倒伏玉米提取结果，共计5种。分别统计5种分类结果中倒伏玉米的栅格数量，然后将栅格数与单个像元的面积作乘求得每种结果中的倒伏玉米面积。

为验证不同特征组合的分类精度，根据目视解译，选取4类典型倒伏作物形态的验证样区各4处，其中未倒伏A类3 599.26 m2、未倒伏B类3 563.09 m2、倒伏3 512.17 m2、阴影区域3 843.3 m2。选取的验证样区与特征统计样区不能有重叠，验证样区的选取结果如图4b，样区选取后目视勾画并统计各类别验证样区中的目标类别面积，以此作为各类别验证的实测面积。

在影像中，倒伏玉米的线性种植结构和垂直生长特征不明显,并且其叶片呈现淡黄色特征[34]，结合目视解译勾画地块内部所有的倒伏玉米边界，解译结果如图4c 所示。经统计边界内的倒伏面积为21 611.91 m2，以此作为实测值，验证本文5种特征组合的倒伏面积提取精度。

图4 样区分布及实测倒伏区域Fig.4 Sample area distribution and measured lodging area

空间一致性是将指定位置的分类结果和验证样本所对应的类别进行比较，多采用混淆矩阵来度量[35]。而Kappa 系数能够揭示景观的空间变化信息，说明数量和位置信息的变换，全面反映分类结果的精度[36]。因此该研究采用Kappa 系数评价倒伏提取结果的空间一致性。计算公式如下：

式中N 是验证样本的总数；n 为混淆矩阵中的总列数；xii为混淆矩阵中第i行、第i列上的样点数；xi+和x+i分别是第i行和第i列的总样点数；K是Kappa系数。

2 结果与分析

2.1 光谱特征差异

计算每类特征统计样区在各波段的反射率均值，绘制4 类作物形态的光谱反射率曲线，如图5 所示，根据图5，与未倒伏A 类相比，倒伏玉米由蓝波段至近红外波段的反射率平均上升0.04；与未倒伏B 类相比，倒伏玉米各波段的反射率平均上升0.01。前者的反射率上升幅度远大于后者，这主要是未倒伏B 类玉米的各波段反射率大于未倒伏A 类玉米的原因。未倒伏A 类与B类之间在绿波段和红边波段上反射率数值差异要高于其他波段。阴影区域在蓝、绿、红波段上与其他3 类区域的反射率差异较大，在红边和近红外波段上反射率差异较小。

图5 4类典型作物形态的光谱反射率曲线Fig.5 Spectral reflectance curves of 4 typical crop morphology

由蓝波段至绿波段，阴影的光谱反射率变化趋于平缓，而其他3 类作物形态的反射率变化趋势相近，且变化速率较快。未倒伏A类与阴影在绿波段至红波段上变化趋势较缓，与未倒伏B类和倒伏形态的反射率变化趋势差异较大。在红波段至近红外波段上，未倒伏A类、未倒伏B类和倒伏形态的光谱反射率变化趋势大体相同，均呈现较快的变化速率。但阴影在红波段至红边波段的反射率变化要大于其他3类作物形态，而在红边波段至近红外波段上与其他3类作物形态的光谱曲线的变化基本相似。

2.2 植被指数差异

计算特征统计样区内各类作物形态的植被指数如表1，在两类植被指数中阴影与其他3 类作物形态之间的差异最大。在NDVI 中未倒伏B 类与倒伏的数值差异最小，而在NDVIR-edge中，差异最小的是未倒伏A 类与未倒伏B 类。对于未倒伏的玉米而言，相机探测的反射率以叶片为主，而玉米在倒伏后，相机所探测的反射率以茎秆和叶片为主。由于茎秆中的叶绿素含量低于叶片，因此玉米发生倒伏后植被指数也会下降，其中NDVI 平均下降了0.08，NDVIR-edge下降0.07。

表1 4类典型作物形态的植被指数Table 1 Vegetation index of 4 typical crop morphology

2.3 纹理特征差异

计算特征统计样区内每类作物形态的纹理特征均值，分析不同形态之间的相对差异度。根据表2，倒伏与2类未倒伏玉米的纹理特征之间有不同的相对差异度。在倒伏与未倒伏A 类玉米的纹理特征对比中，差异最大的是均值纹理，在该特征空间内，2 类作物形态的的相对差异度达1.11，比其他纹理特征平均高0.59；而在倒伏与未倒伏B类玉米的纹理差异分析中，相对差异最大的是对比度，在该特征中，2类作物形态的差异度为2.76，比其他纹理向量平均高0.8；与其他纹理特征相比，二阶矩纹理特征在阴影与倒伏玉米、阴影与未倒伏A类玉米、阴影与未倒伏B类玉米的差异分析中，具有最高的区分度，其平均相对差异为31.87，比其他纹理特征平均高29.63。

表2 纹理特征的相对差异统计Table 2 Relative difference statistics of texture features

2.4 倒伏玉米提取结果与分析

本研究在使用植被指数、多类、单类纹理特征以及光谱反射率特征集进行倒伏提取时，均采用最大似然法。图6a 为利用多类纹理特征集合（均值、对比度以及协同性纹理）提取倒伏玉米的结果；图6b 为使用单类纹理特征集（均值纹理）所得到的倒伏玉米结果；图6c 为利用蓝、绿、红波段的光谱特征识别倒伏玉米的结果；图6d 和6e分别为NDVI和NDVIR-edge特征的倒伏提取结果。

图6 不同特征的玉米倒伏面积提取结果Fig.6 Extraction results of maize lodging area under different features

利用每类作物形态的验证样区，验证不同特征组合的分类误差，统计结果如表3 所示。根据表3 可知，使用多类纹理特征识别阴影和未倒伏A类形态的误差稍高于植被指数和单类纹理特征识别方法，但在对其他作物形态进行分类时，多类纹理特征的分类误差显著低于光谱、植被指数和单类纹理特征识别方法。基于5 种特征组合所提取的倒伏结果如图6 所示，根据此结果计算倒伏面积提取误差，利用100 个倒伏和未倒伏验证样本，计算提取结果的Kappa系数，计算结果如表4所示。

由表4 可知，纹理类的特征组合提取倒伏结果的Kappa 系数最高，平均为0.61，比光谱和植被指数特征结果的Kappa系数分别高0.28、0.46。分析原因是玉米发生倒伏后其主要的变化是形态结构的变化，而光谱和植被指数特征变化相对较小，特别是完熟期未倒伏B 类玉米和倒伏玉米间的光谱和植被指数特征的物理差异很小，因此利用这两类特征进行倒伏提取会产生较多的错分，而纹理特征能够很好地反映未倒伏和倒伏玉米间的空间形态差异，因此利用该类特征进行倒伏提取产生的错分较少，空间一致度较高。

在5 种特征组合中，基于多类纹理特征进行倒伏提取的误差最低，4 类典型作物形态的识别平均误差为9.82%，Kappa系数最高，与单类纹理特征提取结果相比，其Kappa系数提高了0.47，倒伏提取误差降低了27.57 个百分点。分析原因是单类纹理特征虽然能够在一定程度上反映作物的空间域形态信息，但研究区内存在不同倒伏程度和不同生长状况的玉米，它们均具有不同的形态特点，单类纹理的空间映射并不能很好地反映倒伏和未倒伏玉米间的形态差别，而多类纹理特征中所含有的多方向的空间映射信息能够较好地体现倒伏和未倒伏玉米间的形态差异，因此具有较好的空间一致性。

利用光谱特征和NDVI 所提取的倒伏面积过大，其原因主要是大量的未倒伏B 类形态被划分为倒伏，从而使倒伏像元数增多；而利用NDVIR-edge所得的倒伏面积偏大的主要原因是，未倒伏B 类形态和受边缘效应影响的未倒伏玉米被错分为倒伏；基于单类纹理特征进行分类，部分未倒伏B 类形态也被错分为倒伏，但提取的倒伏面积反而过小，分析其原因是更多的倒伏玉米被划分为未倒伏玉米。

表3 不同特征下4类典型作物形态的分类误差Table 3 Classification error of 4 classes of typical crop morphology under different features

表4 不同特征下玉米倒伏面积提取误差Table 4 Extraction error of maize lodging area under different features

3 讨论

本文利用无人机多光谱影像的植被指数、纹理和光谱特征结合最大似然法分类，对完熟期倒伏玉米面积的提取进行了定量分析，得到了以下结果：

1）地块内土壤和作物的时空异质性会使完熟期的玉米处于不同的生长状态，其中倒伏玉米与长势较好的未倒伏玉米的叶片特征和光谱特征很相似，两类玉米的叶片均呈现淡黄色，各波段的反射率平均差异仅为0.01，对这两类作物形态进行分类会产生较多的错分像元，降低倒伏面积提取的精度。单一反射光谱特征或植被指数特征不能准确地区别完熟期倒伏玉米区域。

2）利用高维多类纹理特征能够显著区别地块内不同作物形态之间的差异，增强类内一致性，从而提高倒伏玉米的识别精度。针对完熟期的倒伏玉米识别，该特征组合具有较高的实用价值。

3）阴影区域的判断是完熟期倒伏玉米精准识别的一个难点。未倒伏玉米地块内部的阴影主要包括叶片间隙与防护林阴影两部分，而倒伏发生后地块内部不仅存在以上两部分阴影，还包括未倒伏玉米投影到倒伏玉米上的阴影，这部分阴影集中分布在小面积倒伏区域内部以及片状倒伏边缘，应同样算作倒伏面积。研究发现，除NDVI 和NDVIR-edge外，在选取合适的纹理滤波窗口的条件下，单类、多类纹理特征也能够较好地剔除叶片间隙的阴影，但后者的提取精度更高。对于地块边缘的防护林阴影，本文的5 种特征组合均不能对其进行剔除，考虑可通过实地调查，确定防护林阴影为未倒伏区域之后，根据分类结果建立矢量文件，对阴影区进行手动去除。

4 结 论

本研究利用完熟期玉米的无人机多光谱图像，分析了图像中各类区域的植被、光谱以及纹理特性，针对光谱和纹理特征提取方法分别选取了最优的分类特征，然后利用优选的光谱、植被指数、单类、多类纹理特征向量集结合最大似然分类方法提取了倒伏玉米面积，通过分析分类和提取结果得到以下结论：倒伏玉米与未倒伏A 类玉米的光谱特征差异较大，倒伏玉米与未倒伏B 类玉米的光谱特征差异较小，光谱特征组合无法准确区分倒伏玉米和未倒伏B 类玉米，该特征组合所提取的倒伏玉米面积偏差较大，因此不适于完熟期倒伏玉米面积提取；植被指数特征会使受边缘效应影响的未倒伏玉米错分为倒伏玉米，在该类特征下，倒伏玉米提取结果的空间一致性较差。因此也不适于完熟期倒伏玉米面积提取；对比不同特征分类和倒伏面积提取结果发现，利用均值、对比度及协同性的多类纹理特征组合提取的倒伏玉米面积精度最高，在该类特征下分类的平均误差为9.82%，倒伏面积提取误差为3.40%，Kappa 系数为0.84，比其他特征组合的Kappa 系数平均高0.59，因此多类纹理特征组合可以实现完熟期倒伏玉米面积的准确提取；与光谱特征相比，利用植被指数和纹理特征进行分类能够更好地剔除玉米叶片间隙的阴影，但仍无法准确判断防护林阴影下的玉米是否为倒伏玉米；本研究初步证明了多类纹理特征能够增强倒伏与未倒伏玉米间的空间形态差异，结合不同维数下倒伏提取结果的空间一致性分析能够降低特征信息冗余，防止分类结果的过度拟合，进而提高倒伏提取结果的精度。但受当前无人机续航能力与通讯距离的制约，本研究还无法实现大范围的倒伏面积提取。因此在今后的研究中应尝试利用无人机图像与卫星影像结合的方式来扩大倒伏信息的获取范围。