基于高光谱成像技术的油菜苗期光照胁迫诊断

王怡田，张小敏，姜海益，张延宁，林洋洋，饶秀勤*

（1.浙江大学生物系统工程与食品科学学院，杭州 310058；2.农业农村部农产品产地处理装备重点实验室，杭州 310058；3.浙江大学数学科学学院，杭州 310058）

长江流域是中国油菜主产区，也是世界上规模和开发潜力最大的冬油菜生产带，其油菜常年种植面积和产量均占我国的80%以上[1]。在我国油菜产业生产能力增长缓慢、种植效益偏低的形势下，2017 年，农业部种植业管理司在《全国大宗油料作物生产发展规划（2016—2020年）》[2]中强调，推进长江中下游油菜产区规模化、标准化种植，研究和推广高产栽培技术来提高单产水平，以保障我国油菜产业的可持续发展。

光照在油菜苗的形态建成和生长发育过程中起着关键作用[3−4]，在不同的光照条件下油菜苗形态结构、生化组分含量呈现出不同的变化趋势。在每年10月中下旬适宜移栽期内，长江流域油菜产区易出现寡照天气，影响油菜苗期光合作用的反应过程，导致油菜苗不健壮[5]。研究表明：在弱光条件下油菜叶片光合作用较弱，叶片含水量升高，其生物量积累处于较低水平；而随着光照增强，油菜苗生物量积累增加，可使油菜苗健壮生长[6−7]。因此，研究光照胁迫对油菜苗期生长的影响及光照胁迫诊断方法，对保障我国油菜产业的可持续发展具有重要意义。

作物胁迫诊断通常利用实验室生化方法来进行检测，但该方法存在复杂、耗时、高成本和需专业人员操作等问题，不利于规模化、标准化作业[8]。而机器视觉方法虽然在一定程度上提升了检测速度，简化了操作步骤，但其信息获取具有局限性，准确度不高[9−12]。而高光谱成像技术不仅可以获取待测作物每个像素点的高分辨光谱信息，并且能可视化二维光谱信息，其结合了机器视觉技术和光谱分析技术，能满足更广泛的作物胁迫检测需求[13−14]。当作物受到胁迫时，生化组分含量和内部结构发生改变，产生与正常作物不同的光谱响应[15−16]，因此，研究人员基于高光谱成像技术，通过图谱分析技术捕捉到这种响应差异，进而实现作物胁迫检测。如：RUMPF 等[17]采集了3 种甜菜叶部病害在400～1 050 nm 波长范围内的高光谱图像，提取出9 种光谱指数并建立了支持向量机（support vector machine, SVM）分类器。结果表明，该方法对是否染病的分类准确率达97%，对病害类别的多重分类准确率达86%，对不同病害阶段的分类准确率在65%～90%之间。IORI 等[18]采集了感染诺氏疟原虫的小麦在1 000～1 700 nm 波长范围内的高光谱图像，基于广义最小二乘加权（generalized least squares weighting, GLSW）和主成分分析（principal component analysis,PCA）算法，在接种48 h 后实现了健康和染病组织的检测，且筛选出的波长1 650 nm 可用于感染诺氏疟原虫小麦的早期检测。XIE等[19]采集了番茄叶片早疫病、晚疫病在380～1 023 nm 波长范围内的高光谱图像，基于连续投影算法（successive projection algorithm, SPA）提取特征波长，建立了极限学习机（extreme learning machine,ELM）分类器，其对病害类别的多重分类准确率达97.1%。刘燕德等[20]采集了轻度、中度、重度缺锌的柑橘叶片在400～1 000 nm波长范围内的高光谱图像，基于无信息变量消除（uninformative variable elimination, UVE）算法剔除无关信息，并结合遗传算法（genetic algorithm,GA）和SPA 筛选变量，利用ELM 和最小二乘支持向量机（least squares−support vector machine, LS−SVM）构建柑橘黄龙病判别模型。结果表明，降维特征误判率低于全谱特征，提升了检测效果。

综上所述，基于高光谱成像技术的作物胁迫检测主要集中于全谱分析、特征波长分析等[21−22]，光谱特征的差异不仅反映了胁迫的发生，还提供了胁迫持续时间和严重程度的信息，但胁迫早期检测效果有待进一步提升。本文通过光照胁迫实验分析油菜苗期冠层叶片光谱特征对光照胁迫的响应，对比不同光谱特征提取算法效果，分析光谱特征演化规律，确定合适的光照胁迫敏感特征，从而建立光照胁迫判别模型，以期实现油菜苗期光照胁迫早期诊断。

1 材料与方法

1.1 仪器设备

由于不同波段能反映叶片独特的理化性质，为了全面地获取油菜叶片光谱信息，本实验采用Q285高光谱成像系统（德国Cubert 公司）和SOC710 SWIR高光谱成像系统（美国SOC公司）分别采集油菜叶片在特定波长范围内的高光谱图像。如图1所示，将叶片放置于黑色背景的载物台上，采集时关闭暗箱门，单个样本的总采集时间在1 min 内。其中，Q285 高光谱成像系统的光谱范围为450～998 nm，光谱分辨率为4 nm，镜头焦距为23 mm，空间分辨率为1 000×1 000 像素；SOC710 SWIR 高光谱成像系统的光谱范围为917～1 717 nm，光谱分辨率为2.75 nm，镜头焦距为35 mm，空间分辨率为640×568像素。

图1 油菜叶片高光谱图像采集Fig.1 Hyperspectral image acquisition of rapeseed leaves

1.2 实验设计与样本获取

油菜苗期光照胁迫实验在实验室内进行，油菜品种为浙油50，播种前晒种4 h。将配置的含水量为70%的育苗基质与蛭石按照体积比1∶1均匀混合后装盘，然后将油菜种子播种于72 穴的穴盘中，每穴播种1 粒种子，共播种18 盘，之后放入智能人工气候箱中培养40 d。培养条件：昼温为（25±1）℃；夜温为（15±1）℃；相对湿度为（70±5）%；光周期为光照16 h/黑暗8 h；光密度为350 μmol/（m2·s）。

待幼苗长至两叶一心后，进入光照胁迫实验阶段。根据油菜生产的光密度范围，油菜叶片光饱和点约为600 μmol/（m2·s），光补偿点约为50 μmol/（m2·s）[23]。而在人工发光二极管（LED）光源下，400 μmol/（m2·s）较300 和200 μmol/（m2·s）光密度有利于培育壮苗。光密度较大时，油菜叶片吸收的光能可能超过所能利用的能量，发生光抑制现象；光密度较小时，油菜苗通常呈现高度增加而茎粗减小态势，根系相对不发达，形成的油菜叶片保水能力较差，抵御干旱能力较弱，不利于幼苗生长和成活[24−26]。因此，本实验将油菜幼苗的育苗光密度设置为3 组水平：OD1，210 μmol/（m2·s）；OD2，350 μmol/（m2·s）；OD3，420 μmol/（m2·s）。其中，OD2组为对照组。分别在第1—21天，每天从每组幼苗中摘取20株油菜苗的全部冠层叶片，采用Q285 高光谱成像系统和SOC710 SWIR 高光谱成像系统分别采集在450～998 和917～1 717 nm 范围内的高光谱图像。以每株油菜苗所有冠层叶片的平均光谱作为冠层光谱，构建冠层光谱样本集，样本总数均为1 260 个。其中，OD1、OD2、OD3组各为420个。

1.3 光谱数据预处理

由于高光谱成像系统采集的是整个视场范围内的图像，而我们的感兴趣区域（region of interest,ROI）是冠层叶片区域，因此，需要对叶片和背景进行分割。选取油菜样本在450～998 和917～1 717 nm的高光谱图像中100×100像素大小的冠层叶片区域和背景区域，分别求取平均光谱，得到如图2所示的曲线。当冠层叶片和背景的光谱曲线反射率差异较大时，该波长图像可用作背景分割，最终选择546和1 121 nm 波长图像分别对450～998 和917～1 717 nm高光谱图像进行分割，分割效果如图3所示。

图2 450～998 nm（A）和917～1 717 nm（B）的叶片区域和背景区域平均光谱曲线Fig.2 Mean spectral curves of leaf and background areas at 450-998 nm(A)and 917-1 717 nm(B)

图3 546和1 121 nm波长图像分割结果Fig.3 Segmentation results of 546 and 1 121 nm images

提取冠层叶片的平均光谱，去除两端噪声波段，选择498～950 nm 共114 个波长和939～1 681 nm 共267 个波长进行分析。为了消除高频随机噪声和批次差异，依次进行Savitzky−Golay平滑、标准化和归一化处理。采用主成分−马氏距离法，根据Dixon 准则或Chauvenet 准则剔除异常样本。计算OD1、OD2、OD3各组平均光谱及每个样本光谱与平均光谱之间的马氏距离，按马氏距离从小到大顺序排列；指定检出水平α=0.05，剔除水平α=0.01，异常样本剔除情况如图4 所示。第1 天OD1组第6 个光谱样本和第10天OD3组第3个光谱样本被判别为异常样本并予以剔除。最终有21 d 的样本用于特征提取和油菜苗光照胁迫研究，OD1、OD2、OD3组样本总数分别为419、420、419个。

图4 第1天OD1组（A）和第10天OD3组（B）的马氏距离分布Fig.4 Mahalanobis distance distribution of group OD1 on day 1(A)and group OD3 on day 10(B)

1.4 光照胁迫敏感波段提取

1.4.1 基于光谱反射率显著性分析提取光照胁迫敏感波段

冠层光谱特性由冠层叶片生化组分和结构特性等因素共同决定。光照胁迫下，叶片生化组分含量和结构特性发生改变，会产生区别于正常的光谱响应，最直接表现为各波长光谱反射率的差异。因此，首先基于光谱反射率提取光密度胁迫敏感波长。

按照观测天数将冠层光谱样本集划分为d1～d21族，每族分为OD1、OD2、OD3组。指定显著性水平α=0.05，依次对di族λj波长反射率r（di，λj）进行正态检验和方差齐性检验。若3组均满足正态或轻微偏态分布和方差齐性要求，则采用单因素方差分析，否则采用Mann-Whitney-Wilcoxon（MWW）检验，检验r（di，λj）在3 组之间是否存在显著差异。若单因素方差分析结果存在显著差异，则采用Student-Newman-Keuls（SNK）法进行多重比较，检验3组之间是否两两存在显著差异。若两两存在显著差异，则将λj作为di族光照胁迫敏感波长，统计d1～d21族λj出现的频数。

由于提取的波长频数整体偏小，不同节位叶片光合速率随着生长期推移可能发生不同程度改变，导致光照胁迫敏感波长发生偏移，基于光谱反射率无法提取到全生长期光照胁迫敏感波长。因此，分别提取油菜冠层二叶期至五叶期光照胁迫敏感波长，以保留频数占该生长期所含族数比例不低于75%的波长作为单生长期光照胁迫敏感波长。单生长期光照胁迫敏感波长分布情况如图5所示。基于光谱反射率共提取到1个单生长期光照胁迫敏感波段1 361～1 408 nm以及若干波长。

图5 基于光谱反射率显著性分析提取的单生长期光照胁迫敏感波长分布Fig.5 Distribution of light-stress-sensitive wavelength within single growth period based on the significance analysis of spectral reflectance

1.4.2 基于连续小波变换提取光照胁迫敏感波段

连续小波变换（continuous wavelet transform,CWT）将冠层光谱在不同尺度下转换成一系列的小波系数，表示当前尺度下的小波与所对应波长的相似程度，这种多尺度分解的特性利于更加细致地分析光谱响应差异，故采用连续小波变换提取光照胁迫敏感波段。本文选择Marr 小波母函数，采用1～50连续的小波分解尺度，基于小波系数分析光照胁迫敏感波段。以d1族为例：

1）计算d1族光谱样本1～50分解尺度下的小波系数w（si，λj），得到小波系数矩阵。采用单因素方差分析，检验w（si，λj）在d1族3组之间是否存在显著差异，得到P值矩阵灰度图，P值越小，灰度值越低，差异越显著。如图6 所示，在498～950 和939～1 681 nm波段都存在低灰度P值，即存在小波系数在3组之间差异显著的情况。

图6 498～950 nm（A）和939～1 681 nm（B）d1族的P值矩阵灰度图Fig.6 Gray-level matrix map of P value in the family d1 at 498-950 nm(A)and 939-1 681 nm(B)

2）若存在显著差异，采用SNK 法进行多重比较，检验d1族3组之间是否两两存在显著差异；若两两存在显著差异，则将λj作为d1族光照胁迫敏感波长，将w（si，λj）作为d1族敏感小波系数。如图7 所示，浅蓝色、橙色、黄色区域分别代表3 次比较中有1、2、3 次在0.05 水平上差异显著。图7A 中不存在黄色区域，主要是由于OD2和OD3组差异不显著；图7B中黄色区域主要分布在较高尺度，说明短期内光密度胁迫可能引起某些波段特征发生变化，而不是仅改变了某一波长反射率。

图7 498～950 nm（A）和939～1 681 nm（B）d1族的光照胁迫敏感小波系数图Fig.7 Light-stress-sensitive wavelet coefficient map in the family d1 at 498-950 nm(A)and 939-1 681 nm(B)

类似地，依次对其余20 族进行分析，统计λj出现的频数。以保留频数占总族数比例不低于75%的波长作为全生长期光密度敏感波长，共提取到1个全生长期光照胁迫敏感波段939～978 nm，其分布情况如图8 所示。为减少敏感波长漏检情况，同时提取单生长期光密度敏感波长，结果共提取到2个单生长期光照胁迫敏感波段662～682、984～1 045 nm，以及若干波长，其分布情况如图9所示。

图8 基于CWT提取的全生长期光照胁迫敏感波长分布Fig.8 Distribution of light-stress-sensitive wavelengths in whole growth period based on CWT

图9 498～950 nm（A）和939～1 681 nm（B）波段基于CWT提取的单生长期光照胁迫敏感波长分布Fig.9 Distribution of light-stress-sensitive wavelengths within single growth period based on CWT at 498-950 nm (A) and 939-1 681 nm(B)

综上所述，以上2 种算法共提取到1 个全生长期光照胁迫敏感波段939～978 nm，3个单生长期光照胁迫敏感波段662～682、984～1045、1 361～1 408 nm，以及若干波长。将以上提取的波段和波长用于光谱特征提取算法研究。

1.5 光谱特征提取

1.5.1 样本集划分

如表1 所示，根据实验中培育的油菜苗期生长状况，将冠层光谱样本集划分为C1、C2、C3、C4族，分别对应二叶期至五叶期。

表1 冠层光谱样本集划分Table 1 Canopy spectral sample set division

后续基于1.4节提取到的光照胁迫敏感波段，对划分后的样本集进行特征提取算法研究，比较不同算法在各样本集上所提取的特征波长、小波特征和波段特征的差异性以及区分光照胁迫的能力特征。

1.5.2 特征波长提取算法

对C1～C4族进行SPA 处理，指定特征波长数为1～20，根据均方根误差筛选最优特征波长，具体结果如表2所示。SPA提取的特征波长主要分布在近红外波段处，其中984、1 408 nm 波长多次出现，与水分子O—H键第一泛频、第二泛频吸收带有关。

表2 C1～C4特征波长Table 2 Characteristic wavelengths of C1-C4

1.5.3 小波特征提取算法

采用CWT 算法提取C1～C4族小波特征。为避免特征个数较多和特征间相关性较强，导致矩阵求逆的精度下降、所建判别函数不稳定的问题，采用逐步判别分析（stepwise discriminant analysis,SDA）筛选小波特征，寻找合适的特征子集。基于CWT−SDA算法提取的小波特征结果如表3所示。

表3 C1～C4小波系数Table 3 Wavelet coefficients of C1-C4

2 结果与讨论

2.1 单一特征判别模型

Fisher 分类器是一种结构简单的线性监督分类模型，能让分类结果的类内方差最小而类间散度最大，其对噪声有着良好的鲁棒性。本研究使用Fisher 判别模型，分别对基于SPA 提取的特征波长和CWT−SDA 算法提取的小波特征构建单一特征判别模型SPA−Fisher 和CWT−SDA−Fisher，初步评价特征波长区分在轻度弱光、轻度强光和正常光照下生长的油菜苗的能力。采用无放回随机抽样的5 折交叉验证方法进行SPA−Fisher 和CWT−SDA−Fisher 模型建立和验证，将5 次分类准确率的平均值作为模型总体分类准确率。单次实验中每一组分类准确率的计算公式如下：分类准确率=TP/（TP＋FP）×100%。式中，TP 为该组判别正确的样本数，FP为该组样本被误判成其他组的样本数。

2 种模型的分类准确率如表4 和表5 所示。结果表明：SPA−Fisher 和CWT−SDA−Fisher 模型对油菜各生长期平均总体分类准确率分别为74.47%和77.25%，后者的分类准确率和稳定性略优，但分类效果仍不理想；总体分类准确率随胁迫时间增长和生长期推移，呈现先升后降趋势，在三叶期达到最佳分类效果，说明油菜冠层叶片生长初期和后期对分类准确率的影响较大；2 种模型对OD1组的平均分类准确率都显著高于OD3组，说明低光密度胁迫响应更易被区分。

表4 通过SPA 优选的特征波长在Fisher 判别模型下的光照胁迫判别结果Table 4 Discrimination results of light stress of characteristic wavelengths optimized by SPA under Fisher discriminant model

表5 通过CWT-SDA 优选的小波特征在Fisher 判别模型下的光照胁迫判别结果Table 5 Discrimination results of light stress of wavelet coefficients optimized by CWT-SDA under Fisher discriminant model

2.2 分类模型优化

2.2.1 光谱波段特征提取

上述采用SPA和CWT−SDA算法提取的光谱特征没有明显的波段分布特点，分类效果仍不理想。因此，对全生长期光密度敏感波段939～978 nm 进行特征提取。通过提取的曲线面积和特征角正切值tanθ2项波段特征描述939～978 nm波段光谱反射率和曲线形状差异。

采用3 阶多项式拟合939～978 nm 波段（图10A），并计算各波长处切线斜率（图10B）。以切线斜率最小处970 nm为分界点，将939～978 nm波段划分成2 段。可以看出，在939～970 nm 波段处光谱差异较为明显。

图10 939～978 nm波段3阶多项式拟合曲线（A）和切线斜率变化（B）Fig.10 Fitting curve by third-order polynomial method (A)and changes of tangent slope(B)at 939-978 nm

如图11 所示，分别设置M=939 nm，A=970 nm，首先计算MA的曲线面积，计算公式为

图11 939～978 nm波段特征tan θ示意图Fig.11 Schematic diagram of band characteristic tan θ at 939-978 nm

随后，为避免边界点噪声，M、A各向内取一个点，计算M´、A´处切线，M´、A´处切线交点O构成夹角∠M´OA´，∠M´OA´的正切值tanθ作为939～970 nm波段的第2个特征。

2.2.2 光谱特征演化规律

以3 d 为间隔，将观测期划分为D1—D7，其中D1—D5为育苗期，D6—D7为移栽期。对上述2 个波段特征（曲线面积和tanθ）进行测试，分析光谱特征随时间的演化规律。

1）939～970 nm 特征波段。939～970 nm 曲线面积变化趋势如图12A所示。可以看出：随观测期推移，曲线面积呈现先减小后增大的趋势，该值在移栽期表现为OD3＞OD2＞OD1。该波段反射率大小主要与叶片结构和含水量有关。叶片厚度较大时，反射率较高；含水量较低时，反射率较高。939～970 nm 特征角的正切值tanθ变化趋势如图12B 所示。可以看出：随观测期推移，tanθ呈现先减小后增大的趋势，该值在移栽期表现为OD3＞OD2＞OD1。

图12 939～978 nm波段曲线面积（A）和tan θ（B）随观测期的变化趋势Fig.12 Variation trends of areas under curve(A)and tan θ(B)with observation period at 939-978 nm

2）特征波长。984 nm 反射率变化趋势如图13A 所示。可以看出：在观测期内，984 nm 反射率表现为OD3＞OD2＞OD1。研究表明，984 nm附近为水分子O—H 键第二泛频吸收带，984 nm 反射率与叶片含水量呈负相关。1 408 nm光谱反射率变化趋势如图13B 所示。可以看出：在该波长下光谱反射率变化没有明显的一致趋势，与此波长相关的理化特征在不同光密度下表现出不同的变化趋势。

图13 984 nm（A）和1 408 nm（B）反射率随观测期的变化趋势Fig.13 Variation trends of reflectance with observation period at 984 nm(A)and 1 408 nm(B)

2.3 光照胁迫判别模型

基于2.2 节中的4 个光谱特征，即939～978 nm波段曲线面积、特征角正切值tanθ以及984 和1 408 nm 处反射率，建立多特征融合的Fisher 判别模型，采用无放回随机抽样的5 折交叉验证方法进行模型建立和验证，将5 次分类的平均准确率作为模型总体分类准确率。总体分类准确率和OD1组分类准确率结果如图14所示。结果表明：多特征融合的光照胁迫Fisher 判别模型总体分类准确率为75.00%～95.00%，各生长初期或后期，分类准确率较低，在d20族达到最佳分类效果，准确率为95.00%；该模型对OD1组的分类准确率为80.16%～100.00%，即对弱光胁迫样本检测效果较好。多特征融合的Fisher 判别模型的平均分类准确率为86.88%，显著提升了单类特征光照胁迫检测模型的分类效果。同时，发现单片叶相较于冠层叶而言，干扰因素较少，可能更适于光谱响应分析。

图14 多特征融合的Fisher判别模型总体分类准确率（A）和OD1组分类准确率（B）的变化趋势Fig.14 Variation trends of total classification accuracy (A) and OD1 group classification accuracy (B) by Fisher discriminant model based on multi feature fusion

3 结论

本文以移栽前油菜为研究对象，进行了光密度胁迫实验，研究了基于高光谱成像技术的油菜光照胁迫诊断方法，建立了SPA−Fisher 和CWT−SDA−Fisher 2 种单特征油菜苗光照胁迫判别模型。结果表明，这2 个判别模型对不同生长期油菜的平均分类准确率分别为74.47% 和77.25%，其中CWT−SDA−Fisher 模型的准确率和稳健性更优。进一步筛选出939～978 nm 波段曲线面积、特征角正切值tanθ以及984 和1 408 nm 处反射率4 个特征，建立了多特征融合的Fisher 判别模型，其平均分类准确率为86.88%，显著提升了单类特征光照胁迫检测模型的效果。研究还发现，单片叶相较于冠层叶，干扰因素较少，可能更适于进行光谱响应分析。因此，后续将进一步基于不同真叶研究光密度胁迫检测方法，优化光照胁迫判别模型，从而提高检测效果。