0 引言
【研究意义】据国家统计局统计数据表明,玉米自2012年起成为中国的第一大粮食作物,玉米倒伏是由外界因素引发的玉米茎秆从自然直立状态到永久错位的作物受灾现象,是玉米生产过程中普遍存在的问题之一[1]。玉米倒伏是由品种特性、栽培技术以及土壤、气候等内外因素综合作用的结果[2]。随着农业生产技术的发展和田间管理的改善,玉米倒伏由于其本身的不可控性和难处理性已逐渐成为影响玉米高产、稳产的主要限制因素之一[3]。据统计,由于倒伏导致玉米减产一般可达15%—30%,严重区域甚至达到50%或者绝收[4-7]。玉米倒伏多发生在抽雄期至灌浆期,由于此时处于大风暴雨时节,加上玉米生长旺盛、现代农民过分追求产量而形成的稠密种植和品种等原因,易造成玉米茎秆细、长,纤维素和木质素较少等,而这些都极易导致玉米在大风暴雨的环境下发生倒伏,严重影响玉米籽粒产量、品质和机械收获能力[8]。此外,近年来中国农业保险发展迅速,日益受到农业管理部门、种植户和保险企业的高度重视,推广玉米种植保险的关键在于客观、快速、准确地预测玉米在不同受灾强度下的减产率,有助于及时开展保险理赔和救灾补偿,最大限度地减少种植户受灾损失。因此,及时、客观、快速和定量地监测玉米的倒伏发生范围和受灾程度,定量评估倒伏对玉米产量可能造成的损失,可以为农艺救灾技术扶持、救灾物资优化配置以及农业保险快速理赔提供可靠的信息支撑。【前人研究进展】随着遥感技术的发展,其本身所具有的快速、经济、客观、区域性等优点也为作物倒伏监测、灾情评估方面提供了重要支撑[9-11]。目前国内外不少****在作物倒伏遥感监测方面进行了探索性研究。刘良云等[12]研究发现小麦倒伏后冠层光谱反射率随倒伏角度增加而增加,并通过陆地资源卫星LandSat ETM数据建立归一化植被指数(normalized difference vegetation index,NDVI)与倒伏角度的模型来实现小麦倒伏的遥感监测。由于玉米生长状态受品种、播种时间、播种密度、水肥管理等因素影响较大,其个体和群体的差异一般都要高于小麦、水稻、荞麦等作物,因此,此方法应用于玉米倒伏的遥感监测仍有待考究。杨浩[13]、YANG[14]等先后利用Radarsat-2全极化影像数据,提出了利用雷达极化指数监测小麦倒伏的方法。但雷达数据成像周期较长,成本较高,限制了其在作物倒伏监测方面大面积推广。李宗南等[15-16]先后使用了无人机和Wordview影像进行了小尺度玉米倒伏面积提取。ZHANG等[17]利用地面ASD高光谱数据和小波变换评估了倒伏对玉米品质的影响。MARTINEZ-VAZQUEZ等[18]从植株力学特性、土壤特性和风等与倒伏密切相关的参数出发,模拟小麦倒伏的各种参数。BERRY等[19]通过设置不同的倒伏角度对小麦倒伏产量损失数学模型进行了推算和验证。王猛等[20]验证了倒伏玉米相较于正常玉米的冠层光谱反射率和几种典型植被指数整体上都有不同程度的下降。席吉龙等[21]通过模拟灌浆期玉米不同倒伏程度验证了叶绿素荧光指数（Fv/Fm）和产量随着倒伏程度的加大而降低。【本研究切入点】以上这些研究虽然大多集中在小麦上,但也为玉米倒伏的遥感监测提供了方法上的借鉴。由于玉米发生倒伏后,冠层群体结构产生较大变化,探测视场中的茎叶比及植株组分受光条件的变化直接反映于冠层光谱反射率的变化差异。与此同时,其倒伏前后的长势参量也会发生较大变化,因此,可通过倒伏前后的多时相遥感影像进行长势变化信息监测,进而实现倒伏玉米的识别和灾情严重程度评估。【拟解决的关键问题】利用中分辨率环境减灾小卫星HJ-1B CCD多光谱影像(charge-coupled device)提取多种植被指数,分析倒伏前后的不同倒伏样本的植被指数变化情况,采用相关性分析法筛选最佳倒伏指示性植被指数,构建玉米倒伏灾情遥感监测模型,基于正态（偏正态）统计理论的双阈值划分策略,结合实地调查的玉米倒伏灾情划分倒伏等级,实现区域尺度的玉米倒伏受灾范围和灾情等级遥感制图,为灾情评估和灾后补救提供技术上的支撑。1 材料与方法
1.1 试验区域
本研究于2015年8月在藁城市范围内随机调查进行。藁城市位于河北省石家庄市中东部,地处太行山东麓平原,辖区总面积836 km2,其中耕地约549 km2,是典型的农业耕作区。地理坐标为北纬37°51′18″—38°18′44″,东经114°38′45″—114°58′47″。藁城气候属暖温带半湿润大陆性季风气候。由于地处平原,地势、地貌对气候影响不大,气候因素分布比较均匀,表现为冬冷夏热的气候特点。全境四季分明,春季干燥多风,夏季炎热多雨,秋季气温凉爽,冬季寒冷、雨雪稀少。藁城市夏季粮食作物主要为玉米,同时也是河北省玉米的主产区。1.2 数据获取
2015年8月3日藁城市出现了季节性暴风雨,导致研究区内发生了大面积玉米倒伏。本文选取两期HJ-1B CCD2影像作为本次玉米倒伏遥感监测的主要数据源,其参数如表1,时相分别为7月25日(倒伏前,抽雄期）和8月10日（倒伏后,吐丝期）,且由影像中心经纬度和研究区经纬度计算得到试验区的观测角度约为83°,近似星下点垂直监测,影像质量良好,无云遮盖。据当地农技推广人员分析,玉米此时处于抽雄期至吐丝期,生理生长仍在继续,各种生理指标(如叶面积指数、生物量、株高等)还在上升,而玉米遭受倒伏胁迫后群体结构受到较大破坏,虽然玉米自我修复能力较强,但倒伏与未倒伏之间的玉米群体仍会存在结构和生理指标上的巨大差异,即倒伏与未倒伏之间以及不同倒伏胁迫程度之间的玉米群体会在倒伏前后的长势上存在不同的变化差异,因此,倒伏前后的玉米长势变化差异可用来监测玉米倒伏发生范围及其灾情等级。于8月5日（倒伏后第2天）联合藁城市农技推广中心开展野外调查,在全区范围内随机选择大面积玉米地块进行实地观测,由于HJ-1B CCD2传感器空间分辨率为30 m以及尽量避免混合像元造成的像元误差,样本地块面积不小于33.35 km2,并在样本地块中心位置使用差分GPS记录样本地块的空间坐标,并通过对角线取样法记录中心位置30 m×30 m内的倒伏比例,用于表征倒伏灾情严重程度。其中包括未倒伏地块和倒伏地块,共计47个野外观测样本（其中有一个观测样本在玉米种植区域提取完成后验证精度时,被误分类非玉米区域,故只有46个野外观测样本）,随机选取其中的30个野外实测样本为建模样本,剩余16个为检验样本（建模和检验样本内均包含倒伏和未倒伏地块）。Table 1
表1
表1HJ-1B CCD2影像主要参数
Table 1Main parameters of HJ-1B CCD2 image
| 平台 Platform | 挂载传感器 Mount sensor | 波段 Band | 波段名称 Band name | 光谱范围 Spectral range (μm) | 空间分辨率 Spatial resolution (m) | 幅宽 Width (km) | 重访周期 Revisit cycle (d) |
|---|---|---|---|---|---|---|---|
| HJ-1B | CCD 相机 | 1 | 蓝 Blue | 0.43—0.52 | 30 | 360 | 4 |
| 2 | 绿 Green | 0.52—0.60 | |||||
| 3 | 红 Red | 0.63—0.69 | |||||
| 4 | 近红外 Near infrared | 0.76—0.90 |
新窗口打开
1.3 数据处理与玉米种植区提取方法
两幅HJ-1B CCD2影像处理流程如下：1）辐射定标和波段叠加：利用定标公式中的绝对定标增益系数和偏移量将原始多光谱图像4个波段的DN值(digital number)依次转换为辐亮度,并通过ENVI5.1中的Layer Stacking将定标过的4个波段叠加为一个多波段光谱影像;2）大气校正：采用ENVI5.1中自带的FLAASH模块进行大气校正（由于两期影像时间较近,且查询当天天气可知两天天气状况均为晴朗无云天气,影像质量较高,故采用统一大气参数进行大气校正）,消除大气对地物反射率的影响,得到地表地物的真实反射率值;3）几何精校正：两期HJ星影像均以已知地理信息的TM影像为基准图,使用ENVI5.1中的Image Registration Workflow模块采用3次卷积内插法自动选取控制点50个,并辅用人工修正的方式,对影像进行几何精校正,确保纠正误差小于0.5个像元。玉米种植区提取：由于藁城市地处平原,玉米种植集中且单一,同期混淆作物较少,玉米遥感识别较为简单,在野外样本的支持下,采用最大似然监督分类法对7月25日HJ-1B CCD2影像进行玉米提取,分类精度为97.8%,分类效果较好,满足本文后续研究要求。并使用Arcgis软件做出样图如下,藁城市7月25日HJ-1B CCD2遥感影像（a）和玉米种植空间分布图（b）见图1,其中a图由绿光、近红外、红光波段依次组合而成。
显示原图|下载原图ZIP|生成PPT图17月25日HJ-1B CCD2影像与玉米种植空间分布图
-->Fig. 1HJ-1B CCD2 image on July 25 and spatial distribution of maize in Gaocheng
-->
1.4 分析方法
1.4.1 数理统计方法 对于所获取的玉米倒伏数据首先采用数理统计方法进行数据预处理,主要包括利用Excel 2013进行数据整理和图表制作,运用SPSS 19.0统计分析软件的 Bivariate Correlation（Pearson）模块进行相关性分析和General Linear Model模块进行方差分析与显著性P值计算等。1.4.2 阈值划分方法 采用正态(偏正态)统计理论的双阈值划分策略,提出一种基于倒伏前后比值植被指数变化量的玉米倒伏灾情等级监测方法。首先提取用于建模的野外样本点的倒伏比例和ΔRVI,并计算出倒伏与未倒伏玉米各自的均值、最大值、最小值、标准差等,依据野外实测样本倒伏比例和ΔRVI的关联分析,选取未倒伏的最大ΔRVI值作为划分倒伏与未倒伏玉米的一级阈值Y1。针对倒伏区域,计算倒伏样本地块ΔRVI的均值μ和标准差σ,并划分为[Y1,μ-σ]、[μ-σ,μ]、[μ,μ+σ]、[μ+σ,Max] 4个区间,其中Max为ΔRVI影像里的最大值,这里的μ-σ、μ、μ+σ称为二级阈值Y2,并将这4个区间依次对应于未倒伏、轻度倒伏、中度倒伏、重度倒伏4个倒伏等级,最后将二级阈值划分的未倒伏与一级阈值划分的未倒伏按照统计学的包含关系统一归为未倒伏,技术流程图如图2所示。
显示原图|下载原图ZIP|生成PPT图2玉米倒伏受灾程度遥感监测流程图
-->Fig. 2Flow chart of the the degree of maize affected by lodging monitored via remote sensing monitoring
-->
2 结果
2.1 光谱分析
图3是对预处理之后的8月10号（倒伏后）HJ-1B CCD2多光谱影像经过导入地面GPS点（倒伏和未倒伏玉米采样点）和提取多个影像中可以明显标识的城镇点和水体点,然后在Excel中对每类地物的光谱曲线进行均值化处理并作图。其中城镇建筑物与水体的光谱反射率曲线与玉米的多光谱反射率曲线差异较大,同时也是1.3中的玉米分类精度较高的缘由。针对倒伏和未倒伏玉米,其在可见光波段范围,倒伏玉米要高于未倒伏玉米的反射率。然而在近红外波段（830 nm）,未倒伏玉米略高于倒伏玉米的反射率,主要原因在于玉米倒伏后茎对冠层光谱的贡献比增大,导致倒伏玉米与未倒伏玉米分别在近红外和可见光波段分别出现差异。因此可通过光谱波段组合（即植被指数）进行倒伏程度划分,考虑到植被指数的滞后性,可以通过倒伏前后的植被指数差值来尽量削弱滞后性对倒伏监测的影响[22]。
显示原图|下载原图ZIP|生成PPT图3地物的光谱反射率曲线
-->Fig. 3Spectrum reflectance curve of surface features
-->
2.2 植被指数筛选
通过借鉴前人研究和相关农学知识,本文选取了6种植被指数,分别为：归一化植被指数(normalized difference vegetation index,NDVI）、结构不敏感色素指数（structure insensitive pigment index,SIPI）、增强植被指数（enhanced vegetation index,EVI）、比值植被指数（ratio vegetation index,RVI）、大气阻抗植被指数（atmospherically resistant vegetation index,ARVI）、垂直植被指数（perpendicular vegetation index,PVI）,其具体计算公式如表2。Table 2
表2
表2植被指数计算公式
Table 2Formula of vegetation index
| 植被指数 Vegetation index | 简写 Abbreviation | 计算公式 Formula | 引文 Quotation |
|---|---|---|---|
| 比值植被指数 Ratio vegetation index | RVI | Rnir/Rred | [23] |
| 归一化植被指数 Normalized difference vegetation index | NDVI | (Rnir–Rred)/(Rnir+Rred) | [24] |
| 大气抗阻植被指数 Atmospherically resistant vegetation index | ARVI | (Rnir–(2×Rred–Rblue))/(Rnir+(2×Rred–Rblue)) | [25] |
| 增强型植被指数 Enhanced vegetation index | EVI | 2×(Rnir–Rred)/(Rnir+6×Rred+7.5×Rblue+1) | [26] |
| 垂直植被指数 Perpendicular vegetation index | PVI | (Rnir–b)×cosθ–Rred×sinθ) | [27] |
| 结构不敏感色素指数 Structure insensitive pigment index | SIPI | (Rnir–Rblue)/(Rnir–Rred) | [28] |
新窗口打开
为了解析植被指数变化情况与倒伏灾情程度的响应机理,对HJ-1B CCD2影像进行波段运算得到6种植被指数,采用掩膜处理提取出玉米种植区域内的植被指数;计算并分析倒伏前后的玉米地块的植被指数变化量（ΔVIi,定义为倒伏前玉米植被指数减去倒伏后玉米植被指数所得到的差值）;提取建模野外样本的ΔVIi与其实测玉米倒伏比例进行相关性分析,计算相关系数R,如表3所示。
Table 3
表3
表3植被指数差值与倒伏比例的相关性
Table 3Correlation between vegetation index difference and lodging ratio
| 植被指数 Vegetation index | ΔNDVI | ΔRVI | ΔEVI | ΔSIPI | ΔARVI | ΔPVI |
|---|---|---|---|---|---|---|
| 相关系数 Correlation coefficient (R) | 0.8917** | 0.9377** | 0.2119 | -0.1916 | 0.7652* | -0.3795 |
新窗口打开
由表3可知,比值植被指数(RVI)变化量与玉米倒伏比例的相关性最佳,相关系数可达0.9377,说明ΔRVI具有较好的玉米倒伏受灾程度表征能力;ΔEVI、ΔSIPI、ΔPVI与倒伏比例的相关性均较低,表明这3种植被指数无法准确反演玉米倒伏群体结构的变化情况,主要由于3种植被指数其自身的波段组合造成的,加上此时玉米处于营养生长旺盛期,自我恢复能力较强,所以导致这3种植被指数变化量无法有效区分不同状态的玉米;ARVI是大气抗阻植被指数,通过加入蓝光,重新组合蓝、红、近红外3种波段来消除大气（如气溶胶、雾霾等）的散射、透射,减弱土壤背景的影响,但由于此时的玉米生长旺盛、覆盖度很高,土壤背景对于植被的光谱反射率的影响较小,从而表现为表3中与实测倒伏比例相关性不是很高（R= 0.7652）;ΔNDVI、ΔRVI的相关性均很高,相关系数分别达到了0.8917和0.9377,这是由于两个植被指数均是红、近红外波段的组合,能较好地反映植被覆盖度和作物生长状况的差异,但当作物生物量较高即高植被区域时,NDVI植被指数易受饱和问题限制[29],因此,RVI更适用于植被生长旺盛、具有高覆盖度的植被遥感监测。究其原因在于NDVI是用非线性拉伸的方式增强了近红外（NIR）和红光（R）反射率的对比度,而RVI是通过建立比值模型直接扩大了NIR和R的反射率的对比度,因此当地面为高植被区域时,RVI的相对增加速度要高于NDVI的相对增加速度,即RVI对高植被区的灵敏度要优于NDVI,因此,本文中的ΔRVI与玉米倒伏比例的相关性要高于ΔNDVI与玉米倒伏比例的相关性。综上所述,本文选取RVI作为玉米倒伏遥感监测的敏感植被指数。
2.3 玉米倒伏灾情监测模型构建
针对实地倒伏区域,借鉴田宝明等[30]的倒伏等级划分,本文将玉米实测倒伏比例划为4种倒伏分级。将倒伏比例为91%—100%的玉米划分为重度倒伏;倒伏比例为61%—90%的玉米划分为中度倒伏;倒伏比例为21%—60%的玉米划分为轻度倒伏;倒伏比例为0—20%的玉米划分为未倒伏,这是由于玉米恢复能力较强,有轻微倒伏现象的玉米可以在短时间内得到快速自我恢复。经实地调查发现,研究区内发生倒伏面积较大,不同倒伏状态均存在,并由图4分析可知倒伏比例越大,其ΔRVI值越大。提取用于建模的野外样本点的倒伏比例和ΔRVI,并计算出倒伏与未倒伏玉米各自的均值、最大值、最小值、标准差等,统计数据见表4。对已提取出玉米种植区域后的前后两期HJ星的差值影像在ENVI 5.3中进行统计分析,分析结果见图5,且由图5分析可知其像元整体分布趋势基本符合正态（偏正态）分布。因此,本文采用1.4.2的阈值划分方法,基于正态（偏正态）统计理论的双阈值划分策略,实现对玉米倒伏灾情模型的构建。Table 4
表4
表4野外采样点统计属性表
Table 4Statistical attribute table of the sampling point
| 玉米 Maize | 最小值 Min | 最大值 Max | 均值 Mean | 标准差 Std |
|---|---|---|---|---|
| 未倒伏(No lodging)ΔRVI | 0.0125 | 0.6152 | 0.2243 | 0.2619 |
| 倒伏(Lodging)ΔRVI | 1.0954 | 3.1023 | 1.8986 | 0.5769 |
新窗口打开
显示原图|下载原图ZIP|生成PPT图4ΔRVI与倒伏比例散点图倒伏比例0%为未倒伏,倒伏比例100%为完全倒伏,倒伏比例越大,表示倒伏越严重,ΔRVI定义为倒伏前RVI值减去倒伏后RVI值
-->Fig. 4RVI difference scatter plot with lodging ratioLodging proportion of 0% as no lodging, lodging proportion of 100% as complete lodging, the larger the lodging proportion, the more severe the lodging, ΔRVI is defined as the RVI value prior to lodging minus that of after lodging
-->
显示原图|下载原图ZIP|生成PPT图5ΔRVI影像像元统计图
-->Fig. 5Pixel statistical chart of RVI difference image
-->
2.4 玉米倒伏灾情遥感监测结果
根据玉米倒伏灾情遥感监测方法对ΔRVI影像进行藁城市的玉米倒伏分级,其遥感监测分级图为图6。通过正态（偏正态）统计理论划分的倒伏分级可以清晰地反映出玉米在整个藁城的受灾空间整体分布趋势是藁城中部偏西和西南部倒伏最严重,藁城东南部和藁城最北部倒伏较轻,藁城中部偏北和东部基本未倒伏,该空间分布监测结果与当地农技推广部门掌握的灾情信息基本一致。因此可以根据倒伏分级进行灾情评估或早期灾后扶苗,以减少倒伏对玉米产量的损失。
显示原图|下载原图ZIP|生成PPT图6藁城市玉米倒伏灾情等级遥感监测图
-->Fig. 6Remote sensing monitoring rank chart of maize lodging disaster in Gaocheng
-->
2.5 精度评价
为较真实、客观地实现模型精度的评价,本文采用混淆矩阵的方法进行精度评价。即选取未参与建模的16个野外实测样本点,将其实测倒伏等级与对应GPS位置的倒伏遥感监测等级汇总制成一个混淆矩阵,按照计算公式分别计算其总体分类精度和Kappa系数。为增强精度评价的可靠性,增加试验次数至50次,每次均采用放回随机抽样,每次均计算其总体分类精度和Kappa系数,最后选取50次总体分类精度和Kappa系数的均值作为最终的精度评价指标。经计算得到得其总体分类精度Pc为85.7%,Kappa系数K为0.804,两种精度验证结果均表明本文的倒伏遥感监测模型精度较高,可以实现对区域内不同倒伏状态空间分布的遥感监测。3 讨论
玉米发生倒伏后,受倒伏胁迫的影响,不同倒伏程度的玉米在后期的自我恢复中所需要的时间以及恢复程度均有较大差异,主要表现在倒伏后玉米的不同长势上,而且倒伏灾情越严重,其对后期长势回复的干扰就越大[31],所以可以通过能够表征作物长势的植被指数差值来反演玉米倒伏的受灾程度[14]。本文选取倒伏发生前后的两幅环境减灾小卫星多光谱影像,提取多种植被指数并进行差值运算,通过相关性分析筛选敏感性最高的植被指数差值,基于正态（偏正态）统计理论的双阈值划分策略进行阈值划定,构建基于植被指数变化分析的玉米倒伏灾情遥感监测模型,最终实现藁城市玉米倒伏的遥感监测。此次关于玉米倒伏的研究仅仅只是玉米生长周期内的一个生育期——抽雄期,并不能完全推广整个生育周期内的玉米倒伏研究。由于灌浆期玉米穗粒开始灌浆,营养生长基本停止,加上处于夏季暴风雨的多发时节,所以灌浆期也是玉米极易发生倒伏的生育期[21,32],并且在相同外部条件下玉米发生倒伏的概率以及倒伏后玉米的恢复能力都有了很大的改变[33],为基于长势变化的玉米倒伏遥感监测方法增添了许多不确定性。
通过相关性分析提取的ΔRVI和双阈值划分方法,可以定量地反演玉米不同倒伏程度的空间分布情况,但是由于不同生育期的玉米在生长状态、群体生理指标以及环境参数都存在一定的差异[4],所以仅仅以抽雄期的玉米进行倒伏遥感监测方面的研究,并不能完全的推广到整个玉米生育期的倒伏遥感监测。在今后的研究工作中,可以加入对玉米多个生育期倒伏的遥感监测的研究,以求找到一个最适宜的综合指标和精度更高的阈值选择方法进行整个玉米生育期的遥感监测,以便后期进行灾情评估、灾后补救和产量减损评估。
4 结论
4.1 通过对倒伏前后植被指数变化量与倒伏比例的相关性分析发现,比值植被指数差值相关性最高,达到0.9377,可作为定量反演玉米倒伏遥感监测的最佳植被指数。4.2 玉米发生倒伏后,比值植被指数差值随着倒伏程度的增加而增大。表明倒伏越严重,玉米的长势就越差,从而导致倒伏后比值植被指数变小,比值植被指数差值变大。并且对其差值影像统计分析发现其大体分布上符合正态（偏正态）分布。
4.3 通过提取多种植被指数,可采用基于正态（偏正态）统计理论的双阈值划定策略,构建基于倒伏前后植被指数变化分析的玉米倒伏灾情遥感监测模型,实现玉米倒伏受灾范围以及灾情等级空间分布制图。
4.4 通过构建混淆矩阵,选择总体分类精度和Kappa系数两个指标进行精度验证,总体分类精度为85.7%,Kappa系数为0.804,表明本研究的倒伏遥感监测模型精度较高。
The authors have declared that no competing interests exist.
