基于无人机遥感影像的多品种水稻抽穗时间估算

莫佳才，彭 漪，刘小娟，王 靖

（武汉大学遥感信息工程学院，湖北 武汉 430079）

水稻是我国三大粮食作物之一，近年来我国每年的水稻总产量稳定在2亿t以上，约占全国粮食总产量的30%。水稻产量维系着国家粮食安全，提高水稻总产量仍是我国水稻产业今后的工作重点[1]，而水稻的育种选种是实现这一目标的基础工作之一。传统的水稻育种选种工作依赖于水稻工作者人工观察和记录，费时费力，而遥感技术能够便捷快速地捕捉信息，提高工作效率，在水稻育种栽培上正得到越来越广泛的应用。水稻生育期监测是对水稻各生育时期的生长状态信息进行记载的过程，对水稻的田间管理和产量预估来说至关重要[2]。在移栽稻中，水稻移栽后的生育期可分为返青期、分蘖期、拔节期、孕穗期、抽穗期、乳熟期和成熟期[3]。其中抽穗期是育种者选择优良亲本和进行授粉操作的关键时期，准确地监测抽穗期对于水稻育种和产量估算具有重要意义[4]。

研究表明，在水稻生长发育过程中，其冠层光谱反射率的变化呈现出一定的规律性，由反射率计算得到的遥感植被指数（Vegetation index，VI）的规律性则更为明显[5]。水稻的植被指数通常在营养生长期逐渐上升，到抽穗期达到最高值，随后转为逐渐下降，直至收割[6-8]。正是这种规律性，使研究者们能够从水稻植被指数的时间序列曲线中识别出水稻的生育时期。由于气象条件、观测仪器扰动等内外因素的影响，观测到的原始曲线通常存在噪声，需要进行平滑处理，常用的方法有滤波法、动态阈值法和函数拟合法等[9]。平滑后的时序植被指数曲线存在特征节点，研究者们常从这些节点直接提取水稻抽穗期。如WANG 等[10]在浙江省德清县、杨浩等[11]在湖南省沅江市用国产HJ卫星数据，从时序植被指数最大值识别水稻抽穗期，误差在5 d 以内；MA 等[12]以多个研究区的地面和无人机遥感数据为基础，将拟合的植被指数时序曲线的曲率最大值点对应的时期认为是水稻的抽穗期，取得了满意的效果。在小范围、水稻品种较为单一的情况下，这种方法可以达到一定的精确度，但当水稻品种复杂化时，误差会变大。如SUN 等[13]用MODIS EVI 时序曲线最大值提取全中国不同地区站点的水稻抽穗期，绝大多数结果与实际相差小于16 d；WANG 等[14]和LI 等[15]分别从MODIS NDVI 序列最大值中识别湖南省和江西省水稻抽穗期，均方根误差（RMSE）分别为20 d 与10 d。一般而言，水稻抽穗期时长为10～15 d[16]，因而，上述3 个研究试验提取的抽穗期结果误差过大。相关研究在以下2个方面存在不足。一方面是多基于卫星遥感数据对水稻抽穗期进行识别。与卫星相比，无人机具有方便、灵活和受大气影响较小的特点，使用无人机拍摄遥感影像可以获得更高的空间和时间分辨率，更适用于田间尺度的精细农业，因此需要深入研究基于无人机遥感数据的水稻抽穗期估算方法。另一方面，多从植被指数时序曲线的特征点上直接提取水稻抽穗期，忽略了特征点与水稻抽穗期之间可能存在的数学关系以及水稻品种的多样性。因此，在已有研究的基础上，利用3个独立试验区的1 102个水稻品种的多时相无人机遥感数据，构建多时相植被指数曲线，提取关键特征点，探索其与水稻抽穗期的数学联系，以此实现多品种水稻的抽穗期估算，为多品种水稻育种和大范围农业监测提供参考和指导。

1 材料和方法

1.1 试验区域

本研究试验区域由三部分组成，分别为湖北省鄂州市的试验区1 和试验区2、海南省陵水县的试验区3。鄂州市的2 个试验区（30°22′32″N，114°44′46″E）均位于鄂州市杜山镇东港村，其中试验区1 的试验田分成1 000 余个面积约1.6 m²的小区，共种植1 014个水稻品种，于2018年夏季进行试验；试验区2 的试验田由48 个36 m²的小区组成，种植了48 个水稻品种，于2019 年夏季进行试验。试验区3（18°31′51″N，110°3′18″E）坐落在陵水县椰林镇文官村，包含种植在40 个60 m²小区的40 个水稻品种，于2018 年冬春季进行试验。3 个试验区的位置和田块如图1 所示。3 个试验区之间相互独立且存在明显差异：鄂州市与陵水县分别位于亚热带与热带，气候差异明显，试验期间光温差异大；种植的水稻品种差异大，仅在试验区2 和试验区3 之间存在少量相同品种；水稻品种的抽穗时间差异大。

图1 试验区示意图Fig.1 Location of the experiment

1.2 数据获取与预处理

1.2.1 实际抽穗日期观测 水稻抽穗是水稻在分蘖拔节结束后，完全发育的稻穗逐渐伸出顶部的现象[17]。当稻田中群体抽穗率达到10%时，为水稻的始穗期[18]。在水稻生长发育期间，实地持续观测水稻长势情况，判断并记录水稻所处生育时期，由此获得各试验区每个品种的实际始穗时间。3个试验区的品种、气象和观测情况等信息如表1所示。

表1 各试验区数据概况Tab.1 Data overview of the experiment areas

1.2.2 无人机数据采集与预处理 从水稻返青期至成熟期，每隔一定的时间，使用无人机执行飞行任务，搭载多光谱相机对试验区拍摄遥感影像，获取各小区水稻的遥感数据。试验区1、2、3分别获取了10、13、15 期的无人机数据。本研究使用的是美国Tetracam 公司生产的Micro-MCA 12 通道多光谱相机，内含从可见光到近红外的12 个波段通道，中心波长分别为490、520、550、570、670、680、700、720、800、850、900、950 nm，其中前10个波段的宽度为10 nm，第11、12 个波段的宽度则分别为20、40 nm。为保证影像空间分辨率的一致性与辐射信息的准确性，3 个试验区的无人机采用了基本相同的飞行高度且尽量选择在晴朗无云或少云的天气状况下飞行。对拍摄得到的原始影像进行波段间几何配准，使12 个波段的影像处于同一坐标系下。随后，以影像拍摄期间铺设在水稻田旁的平整路面上的标准定标毯为基准，采用线性模型对影像进行辐射定标，将原始影像各波段的像元值转化为光谱反射率。为获取各品种水稻的反射率，在每个小区内部靠中心位置勾画合适大小的感兴趣区（Region of interest，ROI），计算ROI 内像元的平均反射率，作为该小区水稻光谱反射率的代表值。

1.2.3 植被指数提取 选择了4 种常用植被指数，分别为归一化植被指数（Normalized difference vegetation index，NDVI）、归一化差分红边植被指数（Normalized difference red edge index，NDRE）、简单比值植被指数（Simple ratio index，SR）、二波段增强型植被指数（Two-band enhanced vegetation index，EVI2），对应计算公式见表2。其中，ρNIR、ρred和ρrededge分别表示近红外、红和红边波段的反射率，在本研究中对应采用800、670、720 nm波段。

表2 所用植被指数Tab.2 Vegetation index used

1.3 研究方法

1.3.1 多时相VI 曲线拟合 原始的多时相VI 曲线虽然存在明显的变化规律，但由于无人机数据采集时的气象条件变化、无人机及内部传感器扰动等内外因素的影响，曲线会存在一定的噪音，使多时相曲线不够平滑甚至有突变现象，因此需对原始曲线进行平滑处理。采用函数拟合法来重构多时相VI曲线。常用的拟合函数包括双logistic 函数、非对称高斯函数和对称高斯函数等，其中非对称高斯函数常采用双高斯函数形式，对称高斯函数即单高斯函数。前2 种函数的表达式中，需要解算的未知参数个数均至少有6 个[23]，而本研究试验区的影像期数较少，用这2 种函数易产生过拟合现象。对称高斯函数，即单高斯函数的参数个数最少，为4个。高斯函数的基本公式为：

式中，a是高斯曲线的峰值，b是峰值对应的横坐标，c是标准差，决定高斯钟形曲线的宽度，d为曲线的常数项。考虑到水稻种植前和收割后水稻田中土壤的各波段反射率无明显差别，VI接近于0，故公式中常数项d设置为0，实际为3 个参数，不易过拟合。

1.3.2 拟合曲线特征提取与抽穗时间估算模型构建 水稻抽穗期是水稻从营养生长转为生殖生长的过渡期，植株从只有绿色的叶片发展为绿色叶片和黄绿色稻穗并存，这一变化反映在冠层光谱上。对各试验区全部品种的多时相VI 曲线进行初步探究后发现，实测始穗日期都集中在拟合曲线峰值右侧，实测始穗DAS与曲线峰值DAS之间存在很强的线性相关性。据此，提出本研究的抽穗时间估算模型：

其中，IHDDAS为initial heading date 的DAS，即始穗日期的DAS（简称为始穗DAS），PEAKDAS为峰值的DAS，a1和a2分别为线性函数关系式的斜率和截距。

1.3.3 模型应用与评价 将水稻抽穗时间估算模型首先在建模区应用，随后将模型推广至完全独立的验证区进行检验，并使用Pearson相关系数r、均方根误差（Root mean square error，RMSE）和相对均方根误差（Relative root mean square error，RRMSE）来综合评价此模型的适用性。

2 结果与分析

2.1 VI与多时相曲线拟合方法选择

在计算各时期的植被指数后，可获得各品种的多时相VI序列。图2展示了试验区2部分水稻品种的多时相VI 散点图。由此可见，NDVI 在观测初期就迅速上升到全生育期内的最大值，随后基本保持平缓，到后期略有下降，在中后期的数值很接近，无法以此区分生育期；NDRE 则在整个生育期呈现先上升后下降的明显变化，且曲线基本上左右对称，该曲线形态对于抽穗期的识别十分有利；SR序列与NDVI 相似，且由于水稻封垄后红波段反射率极小，SR 容易出现异常大的值；EVI2 表现优于NDVI 和SR，但波动比NDRE 大。因此，后续研究将选用NDRE继续进行。四者对比可发现由红边波段构建的植被指数的优越性，这与前人的研究结果相符[24]。

图2 试验区2部分水稻品种的多时相VI散点图Fig.2 Scatter diagram of multi-time VIs of some cultivars in experiment area 2

以试验区2 的珞优9311 品种为例，使用3 种拟合函数对该品种的多时相NDRE 曲线进行拟合，结果如图3 所示。原始数据在前期快速上升，在中间时期有小幅波动，但不影响总体趋势为先上升后下降。双高斯函数曲线虽然整体上对原始数据点拟合最贴近，但明显出现双峰现象，与单季水稻的NDRE 曲线只有单峰的事实不相符，这是曲线的过拟合导致的。双logistic 函数曲线形态整体偏向于营养生长期，顶峰与原始数据的顶峰不重合，不能反映真实的NDRE变化情况。单高斯函数曲线虽然贴合的点最少，但与原始数据的曲线形态最为相似，顶峰位置也最为接近，是比另2种函数更合适的选择。因此，本研究将用单高斯函数进行所有品种的多时相NDRE曲线的拟合。

图3 3种曲线拟合方法示例Fig.3 Examples of three curve fitting methods

2.2 模型构建与应用

为探究峰值日期与实测始穗日期之间的函数关系，由于试验区1 的水稻品种最多且实际抽穗时间跨度最大，为保证建立的模型具有通用性，故将其作为建模区。图4 描绘了试验区1 部分品种的NDRE拟合曲线与实测始穗日期的相对位置。对试验区1 的1 014 个水稻品种的多时相NDRE 曲线，用单高斯函数进行拟合，解算参数b，提取峰值DAS，随后与这些品种的实测始穗DAS作比较（图5）。可以发现，实测始穗DAS 与峰值DAS 高度线性相关，两者的Pearson相关系数达到0.906。根据试验区1的1 014 种水稻的实测始穗DAS 与峰值DAS，以最小二乘准则，解算公式（2）的参数a1和a2，得a1=1.3和a2=-18，代入式（2）中，最终可得本研究的水稻抽穗时间估算模型：

图4 部分水稻品种的实测始穗日期与拟合曲线相对位置Fig.4 Relative position of measured initial heading date and fitting curve of some cultivars

图5 建模区峰值DAS与实测始穗DAS散点图Fig.5 Scatter plot of PEAKDAS and measured IHDDAS in the modeling experiment area

将此模型应用到建模区即试验区1 的1 014 个水稻品种，可以得到这些品种的始穗时间估算值。如图6所示，始穗DAS估算值多与实测值接近，二者构成的二维散点都分布在1∶1 线附近。图6 中有个别点偏离1∶1线较远，经统计，始穗DAS估算值与实测值之差在10 d 以上的只有6 个品种，占比为0.6%。绝大多数品种的差值在5 d 以内，品种个数为910，占比达到89.7%。整体上估算值与实测值的RMSE仅为3.15 d，说明抽穗时间估算模型能够很好地估算建模区水稻的始穗时间。

图6 建模区模型估算结果散点图Fig.6 Scatter plot of estimated results by the model in the modeling experiment

2.3 模型推广与验证

为探究水稻抽穗时间估算模型的泛化性能，将模型IHDDAS= 1.3×PEAKDAS- 18 推广至2 个验证区进行验证，其中验证区1 为试验区2，验证区2 为试验区3。如图7所示，该模型在验证区2的应用效果非常好，40 个品种的始穗时间估算RMSE仅为2.99 d，这个误差小于建模区；同时，经统计仅有3个品种的误差大于5 d。而在验证区1，始穗时间估算值小于实测值的现象较为明显，但误差并不很大，48 个品种的RMSE为5.02 d，误差超过10 d 的也仅有1个品种，整体误差可以接受。综合来看，由建模区建立的水稻抽穗时间估算模型在验证区表现良好，能以较小的误差估算2 个验证区内水稻的始穗时间。表3是模型在3个试验区应用的统计结果。

图7 验证区模型估算结果散点图Fig.7 Scatter plot of estimated results by the model in the verification experiment areas

表3 3个试验区的模型应用结果Tab.3 Statistics of model application results in the three experiment areas

3 结论与讨论

3.1 模型适用性讨论

在本研究的3 个试验区中，建模区和验证区进行的试验有很大的差别。从时间和地点上看，建模区和验证区1 处于同一个试验地点，但试验时间不一致，间隔了1 a；建模区和验证区2 的时间、地点均不一致：前者位于华中地区湖北省的亚热带季风气候区，在2018 年夏季进行试验，6—8 月平均气温为25～33 ℃，月均日照长度约为140 h，后者位于华南地区海南省的热带季风气候区，在2018年冬春季进行试验，1—4 月的平均气温为20～27 ℃，月均日照长度约为180 h。从品种上来看，建模区的试验材料有1 014 个水稻品种，涵盖了中国大部分地区典型的杂交水稻品系，但与验证区1 的48 个品种和验证区2的40个品种均不相同，仅验证区1和2中有8个相同品种。

在生育期长短方面，水稻从抽穗到成熟经历的时间约为40 d[25]，因此可用抽穗期的早晚来反映全生育期的长短。建模区水稻的始穗期从播种后59 d 到103 d 不等，其中大部分在80 d 前，生育期较短；验证区2 水稻的始穗期晚于播种后100 d，基本都在110 d 以上，生育期很长；验证区1 水稻的始穗期位于两者之间，生育期也比较长。总的来说，建模区的水稻属于短生育期，验证区的水稻属于长生育期。值得注意的是，水稻的生育期长短与温度、光照等外界环境条件有很大的关系[26]，陵水冬春季较低的温度和长日照在一定程度上延长了水稻的生育期。

建模区的模型在2个与之差异大的验证区都有良好的表现，这表明本研究建立的水稻抽穗时间估算模型对气候、品种、生育期的长短均不敏感，模型具有很好的适用性。未来试验可以在更多不同的时间、地点、品种上应用，以更广泛地探索模型的适用性。

3.2 抽穗时间与多时相VI曲线峰值时间的关系

本研究根据始穗时间与多时相VI 曲线的峰值时间之间的线性函数关系，建立了估算效果良好的水稻抽穗时间估算模型。然而，在已有文献中，最常用的是直接用峰值时间作为抽穗时间的估算值。为探究2 种方法的估算结果差异，用3 个试验区的数据，对比各品种的多时相NDRE 曲线峰值时间与实际始穗时间，发现峰值DAS 整体上普遍小于实测始穗DAS，且实际抽穗时间越晚，这种差异越明显。经计算，用峰值时间直接估算实际抽穗时间的误差远大于本研究提出的抽穗日期估算模型的估算误差。

这一结果与水稻生长发育期间冠层反射率的变化情况密切相关。从分蘖期到拔节孕穗期，水稻植株的叶片不断增多增大，叶片层数的增加使得附加反射率增加，从而引起近红外波段反射率增大，而红边波段反射率没有明显变化，NDRE 在此期间不断增大；进入幼穗分化期，叶片不再增多并需要向穗部供给养分以维持其生长，此时近红外波段反射率开始下降，红边波段依旧变化不大，故NDRE在此期间开始转为减小；到抽穗期，水稻冠层结构因为绿色稻穗从叶鞘中抽出而发生明显变化，对冠层结构变化敏感的红边波段开始上升，近红外波段变化不大，因此抽穗期内NDRE 减小；之后，稻穗逐渐成熟并弯曲，叶片变黄变小，使得红边波段反射率不断增大，近红外波段反射率减小，从而NDRE继续减小[27-28]。因此，从整个生育期来看，NDRE 在幼穗分化阶段已经达到最大值，达到峰值的时间早于抽穗时间，使得峰值时间不能成为抽穗时间的理想估计值，而本研究根据峰值时间与抽穗时间的线性关系建立的模型能更准确地估算水稻的抽穗时间。

本研究提出了一种以无人机遥感为手段的多品种水稻抽穗时间估算方法。以不同气候区域的上千品种水稻的多时相无人机多光谱遥感数据为基础，根据多时相植被指数曲线的峰值时间与实际始穗时间的强线性关系，建立了统一适用的多品种水稻抽穗时间遥感估算模型。模型在建模区1 014种水稻的始穗期估算误差RMSE达到3.15 d，在品种数为48、40 个的2 个验证区的RMSE分别为5.02、2.99 d，整体估算精确度良好。该模型可以不区分气候、地点、品种等差异，为多品种水稻抽穗时间监测提供了有力途径，是育种栽培和田间管理的理想手段。