大田玉米作物系数无人机多光谱遥感估算方法

韩文霆 邵国敏 马代健 ZHANG Huihui 王 毅 牛亚晓

(1.西北农林科技大学机械与电子工程学院，陕西杨凌 712100； 2.西北农林科技大学水土保持研究所，陕西杨凌 712100；3.美国农业部农业研究服务属，柯林斯堡 CO 80526)

0 引言

作物蒸散量(Evapotranspiration，ET)主要由土壤蒸发(Evaporation)和作物蒸腾(Transpiration)两部分组成[1]，是连接生态与水文过程的重要纽带[2]，其快速监测对准确制定和管理大田灌溉制度及提高大田灌溉用水效率有着非常关键的作用。常用的估算蒸散量的方法包括：涡度相关法、蒸渗仪法、土壤水量平衡法、 波文比能量平衡法、Penman-Monteith模型等[3-4]。其中PM模型基于能量平衡、空气动力学原理和冠层阻力来估算ET，具有明确的物理依据，可以清晰地分析ET变化过程及影响机制[2]。但PM模型计算复杂，因此联合国粮农组织(Food and agricultural organization，FAO)提出了参考作物蒸散量(Reference evapotranspiration，ET0)，由ET0和作物系数估算ET[5]。FAO-56作物系数法[6-7]是世界公认的估算作物蒸散量方法，具有操作简便、精度可靠、实用性强等特点，在世界范围内被广泛应用。

目前，国内外研究人员[1, 3, 8-15]利用蒸渗仪数据、涡度相关数据或者水量平衡法获得作物的实际蒸散量ETc，并采用气象数据获得参考蒸散量ET0，结果已证实上述方法具有很好的精度和普适性，其中双作物系数法估算蒸散量的精度更好[4]，冯禹等[9]的黄土高原旱作玉米蒸散估算结果表明，引入叶面积指数计算覆盖度并校正双作物系数的方法能够较好地估算旱作玉米ET(R2=0.871)。

由于作物生长状况、气象条件和水分胁迫等因素的影响，作物系数必须进行校正才能估算实际蒸散量[7]，但及时发现这些外界因素造成的影响并快速调整作物系数曲线较难，并且在时空尺度上，田间作物用水存在差异性，使用推荐的作物系数计算作物蒸散量误差较大。作物系数Kc和作物冠层光谱植被指数均受作物覆盖度、叶面积指数LAI、绿度等因素的影响[16-18]，此关系使遥感技术估算作物系数成为可能。李贺丽等[4]采用地面光谱技术研究了华北平原地区小麦作物系数与植被指数(NDVI、差值植被指数DVI、RVI、SAVI、绿波段叶绿素指数CIgreen、EVI、修正型土壤调整植被指数MASVI和重归一化植被指数RDVI)的关系以及水分和氮素胁迫对其的影响，结果表明不同施氮水平下Kc无明显规律性差异，且植被指数和Kc的相关性较弱。HUNSAKER等[8]采用地面光谱技术研究了美国西南部沙漠地区棉花充分覆盖前和充分覆盖后的基础作物系数Kcb与NDVI和积温GDD的关系。ER-RAKI等[19]基于地面光谱技术研究了采用摩洛哥中心干旱地区小麦的NDVI估算其基础作物系数Kcb和覆盖度fc的方法，结果显示植被指数与Kcb和fc有很高的相关性。GONTIA等[20]基于卫星遥感技术研究了印度西孟加拉邦小麦植被指数(NDVI和SAVI)与其生育期内每月的作物系数Kc的关系。MASELLI等[21]采用意大利中部地区250 m精度的MODIS-NDVI卫星影像估算其植被覆盖度并结合不同下垫面的气象站观测数据估算区域实际蒸散量。

现有研究中，地面光谱技术可以较好地估算基础作物系数Kcb，但难以估算作物系数Kc；卫星遥感技术可以较好地估算区域作物的作物系数，但卫星的拍摄频率和影像分辨率难以满足大田作物的日蒸散量估算要求。无人机遥感技术由于其平台构建容易、运行维护简便、分辨率高、作业周期短等特点，广泛应用于植被信息监测研究中[22-24]，为估算大田作物系数提供了新的解决方案[25]。目前基于无人机遥感技术估算大田作物的作物系数研究较少，且准确估算水分胁迫大田作物的作物系数对制定大田灌溉制度有重要意义。基于此，本文以2017年内蒙古达拉特旗昭君镇实验站大田玉米、土壤、气象等数据为基础，采用经气象因子和作物覆盖度校正后的双作物系数法计算不同水分胁迫与不同生长时期玉米的作物系数，并使用自主研发的无人机多光谱系统航拍玉米的冠层多光谱影像，提取6种植被指数(SR、EVI、VARI、GNDVI、SAVI和NDVI)，研究不同生长时期(快速生长期、生长中期和生长后期)玉米Kc与6种常用植被指数的关系及水分胁迫对其影响，分析无人机多光谱遥感估算玉米Kc的可行性和适用性，为大田玉米作物系数的遥感估算研究提供技术支持。

1 材料与方法

1.1 实验区概况

实验区域位于内蒙古自治区鄂尔多斯市达拉特旗昭君镇(40°26′0.29″N，109°36′25.99″E)，海拔1 010 m，属于典型温带大陆性气候。实验区域玉米(钧凯918)播种时间为2017年5月20日，出苗时间为2017年6月1日，抽穗时间为2017年7月20日，收获时间为2017年9月7日(青储)，生育期(Day after planting，DAP)共110 d。实验数据采集时间为2017年6月26日—8月29日，对应玉米的生长阶段见表1。实验地块(图1中圆形区域)总面积为1.13 hm2，播种深度约5 cm，行距58 cm，株距25 cm，播种时秸秆还田，并按照当地种植习惯施底肥，出苗后根据出苗情况进行补苗，并在玉米拔节期、大喇叭口期分别追加氮肥(240 kg/hm2)。实验阶段总降雨量为44 mm，日平均气温22.9℃，日平均风速0.5 m/s，日平均相对湿度为57.1%，日平均太阳净辐射为145.1 W/m2。生育期内自然降雨难以满足作物需水要求，主要的水分供给方式为喷灌。实验区土壤为砂壤土，其中砂粒(0.05～2 mm)占80.7%、粉粒(0.002～0.05 mm)占13.7%、黏粒(0～0.002 mm)占5.6%，耕层土壤有机质质量比为47.15 g/kg，田间持水率为29%(体积含水率)，土壤容重为1.56 g/cm3。实验地中选取2个扇形区域(图1)，在区域1和区域2内分别选取一块10 m×10 m的方形样地(图1中A和B区域)，并对样地A和样地B的玉米进行长期监测(2017年6月26日—8月29日)。

表1 2017年玉米生长阶段划分Tab.1 Maize growth stage division in 2017

注：玉米生长阶段划分依据参考FAO-56指南[7]。

图1 实验地俯视图(2017年8月15日，无人机RGB影像，5 cm/pixel)Fig.1 Experimental ground view (2017-08-15,unmanned aerial vehicle RGB image, 5 cm/pixel)

1.2 水分胁迫处理

实验地灌溉方式和灌溉量测量方法如图2所示，采用中心支轴喷灌机(图2a)灌溉，在玉米不同生育期进行不同的水分胁迫处理，如表2所示，田间持水量的50%为充分灌溉，采用雨量筒(图2b)测量样地A和样地B的灌溉量和降雨量，每个样地内均匀放置3个雨量筒，取其平均值。

图2 灌溉方式和灌溉量测量方法Fig.2 Irrigation method and irrigation depth measuring method

1.3 无人机系统及影像数据采集方法

采用自主研发的无人机多光谱影像采集系统(图3)，该系统采用开源飞控Pixhawk控制，经纬 M600型机架，最大载量5 kg，最大续航时间30 min。搭载的多光谱相机(图3b)为RedEdge(MicaSense，USA)，该相机有5个波段(表3)，相机焦距为5.5 mm，视场角为47.2°，图像分辨率为1 280像素×960像素。该相机配备了两块3 m×3 m的灰板(GroupVIII，USA)(图3c)及光强传感器(图3d)。光强传感器可校正航拍过程中外界光线的变化对光谱影像造成的影响，灰板具有固定的反射率，如表3所示，可对航拍影像进行反射率校正，生成反射率影像图，并提取植被指数。

表2 实验区实际灌溉量与降雨量Tab.2 Actual irrigation amount and rainfall in experimental area mm

图3 无人机多光谱影像采集系统Fig.3 UAV multispectral images acquisition system

表3 RedEdge多光谱相机参数及灰板对其中心波长的反射率Tab.3 Multispectral camera parameters and reflectivity of gray plate to its center wavelength

选取晴朗无云天气，于11:00—13:00拍摄玉米冠层多光谱影像，拍摄时镜头垂直向下，如图4所示，每隔3 d拍摄一次，飞行高度70 m，地面分辨率为5 cm/pixel，航向和旁向重叠度均为80%，每次航拍采用固定航线。无人机作业时，首先飞至作业高度并拍摄地面两块灰板，然后进行航拍作业。整个生育期共采集20组多光谱影像数据。

图4 无人机航拍作业现场图Fig.4 UAV aerial photography

1.4 地面数据采集方法

(1)气象站数据采集：农业气象站由中国河北清易电子科技有限公司组装和调试，位于实验区西南100 m处(40°25′55.53″N，109°36′22.69″E)，其下垫面为青草，草平面高度保持在12 cm左右。监测参数包括空气温度、相对湿度、风速、太阳净辐射和降雨量，每30 min采集一次数据。

(2)土壤含水率测量：采用干燥法测量实验区玉米的土壤含水率变化。在A和B样地中心点附近随机选取3组样点，测量180 cm土层含水率变化。测量时，按照30 cm间隔对土壤进行分层采样，每3 d采集1次，并在灌溉前、灌溉后和降雨后进行加测，取180 cm厚土层各采样层土壤含水率的平均值。采用线性插值法对土壤含水率进行处理，得到样地A和样地B土壤含水率的逐日数据序列。

(3)玉米生态指标观测：采用随机采样法对玉米的叶面积指数(LAI)、株高等指标进行采集，每7 d测量一次。在样地A和样地B内随机选取10个采样点，采用LAI-2200C型冠层分析仪测量叶面积指数(LAI)，取其平均测量值；在样地A和样地B内随机选取10株具有代表性的玉米，采用米尺测量株高，取其平均测量值。

1.5 作物系数与植被指数计算方法

1.5.1FAO双作物系数法

FAO-56[7]中指出双作物系数法的计算流程。首先将作物的生育期分成4个生长阶段：生长初期、快速生长期、生长中期和生长后期；然后在FAO-56中查找该作物类型在标准情况下(无任何胁迫，田间管理良好，最小相对湿度为45%，冠层上方2 m高度风速为2 m/s)不同生长阶段的基础作物系数推荐值Kcb,tab，可得到标准情况下的基础作物系数逐日曲线；标况下玉米生长初期、生长中期和生长末期的基础作物系数分别取0.15、1.15、0.15。根据实验区气候条件、作物生理参数(作物高度、LAI等)、土壤水分、田间管理等信息对标况下的基础作物系数曲线进行校正，并按照FAO-56给出的算法[7]计算水分胁迫系数Ks和土壤蒸发系数Ke，计算公式为

(1)

Ke=Kc-KsKcb,a

(2)

其中

Kcb,a=Kc,min+(Kcb,adjust-Kc,min)[1-exp (-0.7LAI)]

(3)

(4)

式中Kcb,adjust——不同生长阶段经气象因子校正后的基础作物系数

Kcb,tab——标准情况下的基础作物系数

u2——不同生长阶段冠层上方2 m处的日平均风速，m/s

RH,min——不同生长阶段日平均最小相对湿度，%

Kcb——不同生长阶段经水分胁迫系数校正后的基础作物系数

h——不同生长阶段平均株高，m

Kcb,a——根据作物实际覆盖度校正的基础作物系数

Kc——作物系数

Kc,min——裸土最小作物系数，取值为0.15～0.20

LAI——实测的冠层叶面积指数

1.5.2玉米植被指数提取方法

植被指数可以对地表植被状况进行简单、有效地度量，前人研究中也使用了植被指数建立作物系数的估算模型。常用的植被指数中，NDVI可以估算植被的生理生态参数(覆盖度、叶面积指数、光合有效辐射等)，SR可监测植被覆盖变化且在植被覆盖浓密的情况下效果最好，SAVI可降低土壤背景的影响，EVI可解决植被指数易饱和的问题，VARI可减小光照差异和大气效应的影响，GNDVI对作物的色素变化敏感[26-27]。因为作物系数Kc综合考虑了作物蒸腾和土壤蒸发两方面的影响因素，所以可监测作物系数变化情况的植被指数应能较好地反映植被和土壤的变化情况。

因此本文选取了6种植被指数(NDVI、SR、SAVI、EVI、VARI和GNDVI)，其计算公式见表4。采用Pix4DMapper软件平台对RedEdge相机拍摄的多光谱影像进行几何校正、高斯均值滤波和拼接处理，生成正射影像图，并利用灰板和光照传感器对光谱影像进行反射率校正，生成反射率影像。采用ENVI软件平台裁剪样地A和样地B的反射率影像，并提取上述6种植被指数。计算时，灰板对多光谱相机各波段中心波长的平均反射率分别作为蓝、绿、红、近红外波段的反射率。

表4 植被指数计算公式Tab.4 Vegetation indices

注：RBlue、RGreen、RRed、RNir分别为灰板对RedEdge相机蓝、绿、红、近红外波段的平均反射率。

2 结果与分析

2.1 玉米生理指标

快速生长期至生长后期样地A和样地B玉米株高和LAI如图5和图6所示。

图5 样地A和样地B玉米株高变化Fig.5 Variation of maize height in sample fields A and B

图6 样地A和样地B玉米LAI变化Fig.6 Variation of maize LAI in sample fields A and B

图5中，样地A和样地B在快速生长期玉米高度逐渐增大，且两样地的株高差异较小，生长中期达到最大值，最大高度分别为273 cm和268 cm。图6中，样地A和样地B在快速生长期玉米LAI逐渐增大；生长中期两样地的LAI都先增大再减小，样地A和样地B玉米LAI最大值分别为3.07和2.70，两样地LAI的整体变化规律符合作物的生长规律；生长中期和生长后期由于样地B的灌溉频率和灌溉量小于样地A的值，样地B出现水分胁迫，所以样地B的LAI小于样地A的值。

2.2 植被指数

快速生长期至生长后期样地A和样地B玉米植被指数如图7和图8所示。

图7 样地A植被指数变化曲线Fig.7 Curves of vegetation indices in sample field A

图8 样地B植被指数变化曲线Fig.8 Curves of vegetation indices in sample field B

由图7和图8可知，样地A和样地B在快速生长期植被指数逐渐增大，生长中期和后期植被指数NDVI、EVI、SAVI、GNDVI、VARI有缓慢降低趋势，植被指数SR有较大的减小趋势；生长中期和后期，水分胁迫使得样地B的叶面积指数小于样地A，植被指数SR对作物的覆盖度变化比较敏感，因此样地B的植被指数SR比样地A减小幅度大。

2.3 作物系数

快速生长期至生长后期样地A和样地B玉米的基础作物系数Kcb、土壤蒸发系数Ke、水分胁迫系数Ks和作物系数Kc的逐日变化曲线如图9所示。

图9 快速生长期至生长后期作物系数逐日变化曲线Fig.9 Curves of daily crop coefficient from rapid growth stage to late growth stage

图9a中，样地A和样地B玉米基础作物系数Kcb在快速生长期逐渐增大，生长中期Kcb达到稳定的较大值，生长后期Kcb逐渐降低，该变化规律主要由于作物覆盖度变化所致；生长中期和生长后期，由于水分胁迫的影响，样地B的作物蒸腾受限，所以其Kcb小于样地A的值；生长中期和后期的Kcb变化是土壤水分条件和作物长势状况等因素综合影响的结果。图9b中，快速生长期和生长中期土壤蒸发系数Ke逐渐减小，主要由于作物覆盖度增大限制了土壤蒸发所致；生长后期样地A的Ke有增大趋势，样地B的Ke逐渐减小，主要原因在于生长后期作物覆盖度逐渐降低，充分灌溉条件下有利于土壤蒸发，水分胁迫会限制土壤蒸发；生长后期的Ke变化也是土壤水分条件和作物长势状况等因素综合影响的结果；Ke出现的波动变化主要因为灌溉调整土壤含水率所致。图9c中，快速生长期至生长后期，由于样地A无水分胁迫，所以其水分胁迫系数Ks保持不变，样地B有水分胁迫，随着作物蒸腾和土壤蒸发，其水分胁迫系数逐渐降低。图9d中，样地A的作物系数Kc在快速生长期逐渐增大，生长中期保持相对稳定的较大值，生长后期逐渐减小，此变化规律和FAO-56中标况下的Kc值变化情况基本相符。样地B的Kc在快速生长期有一定的波动变化，是水分胁迫和土壤蒸发综合作用结果；生长中期和生长后期作物系数Kc开始快速下降，主要由于较长时间的水分胁迫(表2)使得水分胁迫系数Ks持续下降所致。上述结果与冯禹等[9]在黄土高原东部地区春玉米试验结果较为一致。

2.4 玉米植被指数与Kc关系模型

快速生长期样地A和样地B玉米植被指数与作物系数Kc的关系如图10、11所示。

图10 快速生长期样地A植被指数与作物系数Kc的关系(n=25)Fig.10 Relationship between vegetation indices and Kc of maize in sample field A in rapid growing stage (n=25)

图11 快速生长期样地B植被指数与作物系数Kc的关系(n=25)Fig.11 Relationship between vegetation indices and Kc of maize in sample field B in rapid growing stage (n=25)

由图10、11可知，快速生长期样地A玉米植被指数与作物系数Kc的相关性较强(R2为0.731 2～0.940 1，p<0.05，RMSE为0.025 0～0.053 0，n=25)，见表5，其中比值植被指数SR与Kc的相关性最好。主要原因在于快速生长期玉米的植被指数(图7)逐渐增大，灌溉充足玉米的Kc也逐渐增大，因此植被指数和作物系数具有较好的相关性。样地B玉米植被指数和作物系数Kc相关性较弱(R2为0.000 2～0.083 0)，因为快速生长期样地B在水分胁迫和作物覆盖度变化的影响下其Ke波动较大，使得Kc波动也较大，而植被指数变化趋势较平滑(图8)，所以两者相关性较弱。

表5 快速生长期样地A玉米植被指数与作物系数Kc的关系Tab.5 Relationship between vegetation indices and Kc of maize in sample field A in rapid growing stage

生长中期至生长后期样地A和样地B玉米植被指数与作物系数Kc的关系如图12、13所示。

图12 生长中期至生长后期样地A植被指数与作物系数Kc的关系(n=40)Fig.12 Relationship between vegetation indices and Kc of maize in sample field A from middle growing stage to later growing stage (n=40)

图13 生长中期至生长后期样地B植被指数与作物系数Kc的关系(n=40)Fig.13 Relationship between vegetation indices and Kc of maize in sample field B from middle growing stage to later growing stage (n=40)

由图12、13可知，生长中期至生长后期样地A玉米植被指数与作物系数Kc相关性(R2为0.276 5～0.373 2)较弱，因为生长中期至生长后期样地A玉米的覆盖度较高，植被指数已达到饱和，对作物系数变化不敏感，所以两者的相关性较弱。样地B部分植被指数与作物系数Kc的相关性(R2为0.366 2～0.848 7，p<0.05，RMSE为0.114 2～0.233 8，n=40)较强，见表6，主要因为该生长中期至生长后期样地B有水分胁迫，玉米覆盖度逐渐减小，部分植被指数减小趋势较明显，所以两者相关性较高。

表6 生长中期至生长后期样地B植被指数与作物系数Kc的关系Tab.6 Relationship between vegetation indices and crop coefficient Kc of maize in sample field B from middle growing stage to later growing stage

由以上可知，不同生长时期玉米的比值植被指数(SR)和作物系数相关性最好(R2最大为0.940 1，p<0.05，RMSE=0.025 0，n=25)，本研究相对于前人研究取得了较好的结果，李贺丽等[4]基于地面光谱技术建立小麦的VIs-Kc模型，决定系数R2最大为0.15，p<0.01，n=195，PARK等[34]基于卫星遥感技术建立NDVI、LAI和土壤含水率与Kc模型，决定系数R2最大为0.64。本文采用的水分胁迫程度较大，模型的适用性有待进一步验证，后期实验应划分多个水分胁迫水平，对多种水分胁迫玉米的Kc估算精度进行评价。时间上非连续采集土壤含水率，对模型也存在一定的影响。目前水分胁迫的监测手段有热红外影像[35]、地面传感器网络等。后期实验可同时采集多光谱影像和热红外影像数据，并搭建地面传感器网络连续监测作物、土壤等相关参数，对大田玉米作物系数的无人机遥感技术估算模型进一步修正。

3 结论

(1)采用无人机多光谱技术估算Kc具有一定的可行性。在快速生长期充分灌溉和生长后期水分胁迫条件下，均可以采用无人机多光谱技术估算玉米作物系数Kc。

(2)不同时期玉米的植被指数和作物系数的相关性不同。快速生长期充分灌溉玉米的VIs-Kc模型的相关性(R2为0.731 2～0.940 1，p<0.05，RMSE为0.025 0～0.053 0，n=25)与生长中期至生长后期充分灌溉玉米的VIs-Kc模型的相关性(R2为0.276 5～0.373 2，p<0.05，RMSE为0.073 8～0.079 3，n=40)不同；快速生长期水分胁迫玉米的VIs-Kc模型的相关性(R2为0.000 2～0.083 0，p<0.05，RMSE为0.099 4～0.114 4，n=25)与生长中期至生长后期水分胁迫玉米VIs-Kc模型的相关性(R2为0.366 2～0.848 7，p<0.05，RMSE为0.114 2～0.233 8，n=40)不同。

(3)不同水分胁迫玉米的植被指数和作物系数的相关性不同。快速生长期充分灌溉玉米的VIs-Kc模型相关性(R2为0.731 2～0.940 1，p<0.05，n=25)比水分胁迫玉米的VIs-Kc模型的相关性(R2为0.000 2～0.083 0，p<0.05，n=25)强；生长中期至生长后期充分灌溉玉米的VIs-Kc模型相关性(R2为0.276 5～0.373 2，p<0.05，n=40)比水分胁迫玉米的VIs-Kc模型相关性(R2为0.366 2～0.848 7，p<0.05，n=40)弱。

(4)不同植被指数和作物系数相关性不同，比值植被指数SR与作物系数Kc的相关性最好。快速生长期充分灌溉玉米的VIs-Kc模型的相关性由大到小依次为：SR、EVI、VARI、GNDVI、SAVI、NDVI；生长中期至生长后期水分胁迫玉米的VIs-Kc模型的相关性由大到小依次为：SR、GNDVI、VARI、NDVI、SAVI、EVI；其中比值植被指数SR与作物系数Kc的相关性最好。