基于分数阶微分的葡萄叶片SPAD值高光谱遥感反演研究

郭 松,舒 田,刘春艳,冯恩英,王文静,蒋丹垚

(1.贵州省农业科技信息研究所,贵阳 550006;2.西北农林科技大学资源环境学院,陕西 杨凌 712100)

【研究意义】葡萄是世界四大水果之一,其栽培面积和产量居世界首位。葡萄果实不仅是富含氨基酸和矿物质等多种营养物质的鲜食水果,还是酿酒的重要原料[1]。我国葡萄种植区分布广泛,主要有西南山地区、关中平原区和东北高原区等。作为乡村振兴的重要构成部分,葡萄产业已成为促进我国农业经济发展和优化农业种植结构的主要水果产业之一[2],随着葡萄产量需求的不断提升,种植面积的不断扩大,如何快速监测葡萄长势状态已成为当前有待解决的重点问题。传统的人工实地监测不仅费时费力、成本高昂,且效率低下、信息滞后,难以为田间高效管理提供数据支持。近年来,高光谱遥感技术以数据量丰富且数据维度广的特点逐渐应用于植被长势的监测研究[3],其基本原理是定性或定量化描述由于植被长势不同所造成的样本高光谱曲线反射率差异,并将该描述模型化以扩展到未知样本上,从而实现无损且非接触式植被长势的大面积监测。开展基于分数阶微分的葡萄叶片SPAD值高光谱遥感反演研究对我国葡萄种植业信息化管理具有重要意义。【前人研究进展】有关利用高光谱遥感技术监测葡萄长势方面的研究已有文献报道。张旭等[4]以卷积平滑、多元散射校正及一阶导数变换等方法处理葡萄叶片和果实的非成像高光谱数据,使用偏最小二乘法建立葡萄可溶性固体含量估算模型发现,基于果实光谱的最优估算模型优于叶片光谱,其建模与验证的R2分别达0.93和0.86以上。王浩淼等[5]通过无人机搭载高光谱相机获取葡萄园冠层高光谱影像并计算相应影像的葡萄园NDVI空间分布,以NDVI值大于0.75、连续性小于0.60分别作为葡萄长势和缺苗的评价标准,结果表明,葡萄冠层无人机高光谱影像在监测葡萄长势和缺株定位方面具有较高潜力。汤森林等[6]提出一种基于特征选择的长短时记忆神经网络模型用于估算葡萄冠层叶面积指数,该模型精度优于传统偏最小二乘模型和支持向量回归模型,其均方根误差低至0.0745。Wei等[7]对葡萄叶片高光谱曲线采用多种数学变换方式,从不同角度构建多个葡萄茎杆水势估测模型,其中变换光谱结合偏最小二乘构建的模型精度最佳,建模RMSE低至110 kPa。Ryckewaert等[8]利用盆栽葡萄的可见光-近红外高光谱影像估测叶片蒸腾系数和气孔导度,建立序列正交偏最小二乘模型发现,偏最小二乘模型对不同品种葡萄叶片蒸腾系数和气孔导度的预测性能均良好,模型R2在0.6以上;Lyu等[9]采用支持向量回归、随机森林回归以及偏最小二乘法构建葡萄叶片氮、磷、钾、钙和镁含量高光谱估测模型,其中偏最小二乘法表现最佳,钾估测模型精度最高(R2为0.70),钙和镁模型精度较低(R2分别为0.15和0.43)。【本研究切入点】当下基于高光谱技术的葡萄长势研究已较为全面,监测指标包括叶片叶面积指数、氮磷钾含量和水分含量等生理生化参数,模型构建方法涉及传统回归和机器学习,光谱数据获取亦由地面向无人机扩展,但仍存在可改进之处:①光谱变换作为消除背景噪声有效手段之一,当前研究多聚焦于整数阶微分变换和普通数学变换等,而有关分数阶微分变换的研究鲜见报道。与整数阶微分变换相比,分数阶微分变换在消除背景噪声的同时能更大程度地保留光谱有效信息[10-11],同时相较于其他传统光谱变换方式,分数阶微分变换因其步长设置的多样性而具备更为丰富的光谱衍生数据[12]。②基于排列组合的新光谱指数在葡萄长势参数反演方面的研究较为缺乏,基于排列组合的新光谱指数可充分发挥高光谱波段数多而窄的优势,且构建的光谱指数与待反演指标的相关性比传统植被指数的更好[13]。【拟解决的关键问题】以贵州息烽葡萄种植基地的葡萄作为研究对象,结合葡萄冠层叶片高光谱数据和田间实测叶片叶绿素含量(SPAD值),构建基于特征波段和光谱指数的葡萄叶片叶绿素含量反演单因素模型和基于连续投影算法的多因素模型,探究最佳分数阶变换及最优葡萄叶片叶绿素含量估算模型构建方法,为当地及类似地区无损、快速且大面积监测葡萄长势提供参考。

1 材料与方法

1.1 研究区概况

研究区位于贵州省息烽县小寨坝镇红岩村(106°4′1″E,27°8′28″N)(图1),地处黔中山原丘陵中部,平均海拔1000 m,是典型喀斯特地貌区,属亚热带季风气候,冬季温和少雨,夏季高温多雨,年均气温15 ℃,年均降水量1100 mm,集中于5—9月,全年无霜期>260 d,土壤类型为泥质黄壤,适宜果树种植。试验开展于2023年7月25日,葡萄处于成熟期。

图1 研究区地理位置Fig.1 Geographical location of the study area

1.2 试验材料

1.2.1 葡萄叶片 采样对象位于息烽县小寨坝镇红岩村葡萄种植基地,品种为水晶葡萄,树龄15年。

1.2.2 仪器设备 ASD Field Spec3野外便携式地物光谱仪,美国ASD公司;SPAD(Soil plant analysis development)手持式叶绿素仪,日本KONICA MINOLTA公司。

1.3 试验方法

1.3.1 指标测定 (1)叶片高光谱数据。通过ASD Field Spec3野外便携式地物光谱仪采集叶片高光谱数据,该仪器测定波长为350～2500 nm,光谱采样间隔在波长350～1000 nm时为1.4 nm,波长1000～2500 nm时为2.0 nm,选择11:00—14:00(天气晴朗无风、太阳辐射充沛时)进行采样,采样时探头垂直于叶片[14],二者相距30～50 cm。为确保数据准确性,每隔15 min校正1次黑白板;研究区共设置74个采样点:首先根据研究区内各葡萄种植地块面积大小占总面积的百分比加权分配采样点个数,然后在各地块中选择生长状态较好的葡萄植株作为采样点,同一地块不同采样点间遵循至少相隔3 m的原则,每个采样点按照长势梯度选择无病害的嫩叶、成熟叶及老叶各1片进行测量,每片叶子采集10条光谱曲线,按时间序列剔除首尾2条光谱曲线,剩余8条光谱曲线通过Savitzky-Golay二阶平滑滤波后取平均作为目标叶片代表光谱曲线。研究区共采集222条光谱曲线。

(2)叶片叶绿素含量。葡萄叶片叶绿素含量由SPAD手持式叶绿素仪测定[15-16]。SPAD值的测定与光谱数据采集同步进行,分别在每片叶子的尖部和中部测定2个SPAD值,每片叶子获取4个SPAD值,取平均作为目标叶片代表SPAD值。研究区共采集222个SPAD值。

1.3.2 数据预处理 由于植被叶片的光谱响应波段集中于可见光-近红外,故将叶片光谱曲线按1 nm间隔在400～1000 nm重采样[17]。采用Excel 2016对获取的222个葡萄叶片SPAD值以升序排序,按3∶1比例分层抽样,其中167个作为建模集,55个作为验证集(表1)。与总体相比,建模集和验证集各项统计特征差别较小,说明样本划分较为合理,总体、建模集及验证集变异系数均大于30%,属于强变异,原因在于测定SPAD值的葡萄叶片长势差别明显,SPAD值跨度较大[18]。

表1 不同样本集的统计特征Table 1 The statistical characteristics of different sample sets

1.3.3 模型建立 在Matlab 2016a环境下以0.2阶为步长,对成熟期葡萄叶片光谱曲线进行0～1.4阶分数阶微分变换,定义不同微分变换下与SPAD值相关性最高的波段为特征波段,计算不同微分变换下两两波段组合的差值光谱指数(Difference spectral index,DSI)、归一化光谱指数(Normalized difference spectral Index,NDSI)及比值光谱指数(Ratio spectral index,RSI)作为新光谱指数[19],将原始光谱和分数阶微分光谱新光谱指数与相应特征波段作为建立单因素模型的建模变量,选择不同分数阶下相关性最高的单因素模型建模变量,建立基于一元线性函数、幂函数、指数函数以及对数函数的单因素模型。利用连续投影算法(Successive projections algorithm,SPA)不同分数阶微分变换下光谱特征作为自变量,构建基于多元线性回归(Multiple linear regression,MLR)和XGBoost回归(eXtreme gradient boosting regression,XGBR)的多因素模型。其中连续投影算法是一种前向特征选择算法,可同时兼顾自变量对响应变量的解释性以及自变量间的独立性,不仅能有效避免维度灾难(Curse of dimensionality),还能明确自变量物理意义[20]。其基本原理为:定义光谱序列X、响应变量Y以及目标波段数量n,其中X包含m个波段。第一步从X中任意提取1个波段Xi,计算Xi相对其余波段的投影并提取最大投影的波段作为入选波长;第二步更换Xi并重复第一步,直至选入波段数量等于n。由于第一步具有随机性,故每次选定自变量序列后计算一次选定波段对响应变量的预测能力,通常以均方根误差(RMSE)作为评价标准,选择RMSE最小的自变量组合作为最终降维结果。将相关系数绝对值大于0.5的波段以及葡萄叶片SPAD值分别作为待提取自变量集和响应变量,通过SPA不断循环迭代,得到多因素模型建模变量。XGBoost是对GBDT(Gradient boosting decision tree)算法的继承和发展,其基本思路是“集思广益”,即整合众多弱学习器以形成强学习器,弱学习器被分散到不同的决策树中,每添加1棵决策树,就增加1个弱学习器(函数),增加决策树的同时也会学习到目标更多的特征。当模型训练完成后,对象特征会被每棵树分担,待预测目标值即为相应树的反馈值总和[21]。由郑智康等[22-24]的研究结果可知,机器学习算法在拟合非线性关系时优于传统回归,故本研究以不同类型光谱下多因素模型建模变量和葡萄叶片SPAD值为对象,采用XGBoost算法构建SPAD值预测模型,拟合函数选择“reg:linear”。为节约调参时间成本,拟合过程采用遗传算法(Genetic algorithm,GA)寻优,该算法原理是将参数寻优问题转换为自然界中生物基因交叉变异问题,通过模拟自然过程以探索最优解,使得XGBoost拟合过程向着nRMSE尽量小而R2尽量大的方向前进。采用决定系数(R2)和归一化均方根误差(nRMSE)作为模型精度评价指标。二者值域分别为[0,1]和[0,+∞),其中,R2越接近1,表示自变量对响应变量的解释能力越强,模型精度越高;nRMSE越小,表示模型预测值与实测值差异越小,模型预测能力越好[25]。

RSI=Rm/Rn

DSI=Rm-Rn

NDSI=(Rm-Rn)/(Rm+Rn)

2 结果与分析

2.1 葡萄叶片高光谱特征

2.1.1 原始光谱 将222个葡萄叶片光谱曲线样本以SPAD值升序为原则均匀分为3组:各组SPAD值最大值、最小值及平均值分别为(34.09,11.52,23.42)、(46.45,34.32,40.78)和(54.45,46.59,49.59),对每组SPAD值及相应光谱反射率取平均,得到不同SPAD值代表光谱曲线(图2)。总体上,不同SPAD值葡萄叶片原始光谱曲线特征基本一致,形似根号,在可见光区域(400～680 nm),绿色植被吸收蓝光和红光而反射绿光的特性使得光谱曲线上形成蓝谷、红谷及绿峰;在近红外区域(760～1000 nm),叶片细胞内部栅栏组织对近红外光的强烈反射导致高反射平台形成,而红谷与高反射平台之间通过红边(680～760 nm)衔接。局部上,不同SPAD值葡萄叶片原始光谱曲线特征有所不同,在可见光区域,随着SPAD值的增大,整体反射率变小,同时红边向长波方向移动,表征叶片对可见光的吸收整体增强;在近红外区域,随着SPAD值的增大,叶片整体反射率增加。说明,SPAD值越大,葡萄叶片细胞内部栅栏组织越稳定,对近红外光的反射更剧烈。

图2 葡萄叶片的原始光谱特征Fig.2 The original spectral characteristics of grape leaves

2.1.2 分数阶微分光谱 光谱微分变换能突出显示原始光谱曲线陡峭的位置而淡化平缓的位置[24]。从分数阶(图3)看,随着阶数上升,分数阶变换光谱与原始光谱差异越来越大,整体上微分值逐渐向0值靠近,局部上在可见光区域逐渐折线化,在近红外区域则整体下移,其中可见光区域微分阈值由0.2阶下的0.58～5.59逐渐减小至1.4阶下的-0.26～0.11,近红外区域微分阈值则从0.2阶处的12.72～19.60减小至1.4阶处的-0.074～0.028,以0.6阶为分界点,大于0.6阶后各微分曲线开始出现负值;从SPAD值看,随着SPAD值的减小,不同SPAD值分数阶微分光谱曲线在红边706 nm处的极值相对大小由0.2阶下的SPAD值49.59最大、SPAD值40.78次之、SPAD值23.42最小逐步变更为1.4阶下的SPAD值23.42最大、SPAD值40.78次之、SPAD值23.42最小,其中,0.2阶下706 nm处不同SPAD值曲线微分值分别为9.92、8.28和7.42,相应的1.4阶下为0.23、0.27和0.28;同时,当阶数大于0.8阶时,不同SPAD值的微分曲线在可见光区域和部分近红外区域(950～1000 nm)的一致性降低,随着阶数的上升,噪声现象愈加严重。

A、B、C、D、E、F、G分别表示0.2、0.4、0.6、0.8、1.0、1.2和1.4阶下的光谱。A, B, C, D, E, F and G represent the spectrum of order 0.2, 0.4, 0.6, 0.8, 1.0, 1.2 and 1.4, respectively.图3 葡萄叶片的分数阶微分光谱特征Fig.3 Fractional differential spectral characteristics of grape leaves

2.2 光谱指数、特征波段与SPAD值的相关性

根据构建的各类型光谱下单因素模型建模变量及其相关系数(表2),通过查阅相关系数表[26]可知,自由度为221时,所有阶数下单因素建模变量与SPAD值均达极显著相关,相关系数为0.615～0.785;从相关系数看,不同分数阶光谱下的最优新光谱指数均为DSI,随着阶数上升,特征波段、DSI与SPAD值的相关系数呈先升后降趋势,二者相关系数分别在0.6和0.4阶达最大,为0.777和0.785;除0.6与0.8阶外,其余分数阶微分光谱均为特征波段相关性优于DSI相关性。从波段位置看,当阶数递增,特征波段由近红外向红边靠近,而DSI计算波段则由近红外与绿光组合变更为近红外与红边组合,所有单因素模型建模变量中,特征波段与DSI相关系数最高的分别是0.6阶下R(959)和0.4阶下DSI(992,439),二者相应的相关系数分别为0.777和0.785。选择相关性最高的光谱指数形成热图(图4),由此可见,随着阶数上升,DSI相关系数较高的区域逐渐由连片分布转变为零星分布,其中0.0～0.6阶该区域主要集中分布于x、y∈[740 nm,1000 nm],0.8～1.4阶则零星分布于x=700 nm以及y=700 nm两端。

表2 葡萄叶片SPAD值单因素建模的变量Table 2 Single-factor modeling variables of SPAD value in grape leaves

A、B、C、D、E、F、G、H分别表示0.0、0.2、0.4、0.6、0.8、1.0、1.2和1.4阶下的DSI热图。A, B, C, D, E, F,G and H represent the DSI heat map of order 0.0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.2 and 1.4, respectively.图4 葡萄叶片光谱指数与SPAD值的相关性热图Fig.4 Heatmap of correlation between spectral index and SPAD value of grape leaves

2.3 葡萄叶片SPAD值的单因素反演模型

总体上单因素模型包括线性与非线性2种(表3),所有模型建模与验证的R2为0.45～0.65,归一化均方根误差(nRMSE)为19.03%～23.03%。随着分数阶增大,以0.2阶为界,模型精度呈先升后降趋势,单因素模型中以0.4阶下基于DSI的幂函数模型精度最佳,其建模和验证的R2分别为0.65和0.61,相应的nRMSE为19.03%和20.60%,当使用单因素模型反演葡萄叶片SPAD值时应优先考虑此模型。

表3 葡萄叶片SPAD值的单因素反演模型Table 3 Single-factor inversion models of SPAD values of grape leaves

2.4 葡萄叶片SPAD值的多因素反演模型

2.4.1 多因素模型建模变量获取 从表4看出,各分数阶下基于连续投影算法提取的建模变量在5～16个,1.4阶最少而0.6阶最多;选择比率为1.41%～16.09%,降维效果明显。从不同阶数看,除1.4阶外,其余分数阶微分光谱建模变量个数均多于原始光谱,说明分数阶微分变换能减弱不同波段反射率的共线性;从建模变量看,0.0阶主要包含绿光、红光与近红外区域,0.2～0.6阶在0.0阶基础上扩展了绿光区域,相应的0.8～1.4阶则渐渐剔除绿光与近红外区域,说明随着阶数的上升,各波段有效信息逐步向绿光、红光以及红边波段靠拢。

表4 葡萄叶片SPAD值多因素模型建模变量Table 4 Modeling variables of multi-factor model for SPAD of grape leaves

2.4.2 葡萄叶片SPAD值的SPA-MLR模型与其精度 由于SPA选择多因素模型建模变量时已规避了自变量间的共线性,故建立多元线性回归(MLR)模型时不再进行多变量逐步回归。建模结果(表5)表明,与单因素模型相比,各模型建模与验证的R2分别提高0.01～0.15和0.05～0.19,相应nRMSE则分别降低0.61%～3.92%和1.52%～4.79%,可见SPA-MLR模型精度更佳,其中最优模型位于0.6阶微分光谱,其建模与验证的R2分别高达0.71和0.68,nRMSE则分别低至16.91%和17.37%。

表5 葡萄叶片SPAD值的多元线性回归模型Table 5 Multiple linear regression models of SPAD values of grape leaves

2.4.3 葡萄叶片SPAD值的SPA-GA-XGBoost模型与其精度 GA-XGBoost算法拟合结果(表6)表明,与SPA-MLR模型相比,不同分数阶SPA-GA-XGBoost回归模型精度均有所提高,其建模和验证的R2分别提高0.03～0.09和0.03～0.08,相应的nRMSE则降低1.14%～2.46%和0.60%～1.98%;精度最高的SPA-GA-XGBoost模型位于0.6阶,建模和验证的精度分别为(0.79,14.45%)和(0.75,15.54%)。

表6 葡萄叶片SPAD值的SPA-GA-XGBoost回归模型Table 6 SPA-GA-XGBoost regression models of SPAD values of grape leaves

3 讨 论

定性化描述是反演模型定量化表达的前提,不同SPAD值葡萄叶片原始光谱总体一致而局部不同,光谱曲线细节上的差异在数理统计方面体现为反射率与SPAD值的正、负相关关系,随着SPAD值的上升,反射率在可见光区域整体下降而在近红外区域整体上升,表征葡萄叶片光合色素含量增加且细胞栅栏组织更为稳定,该结论与王婷婷[27]的研究结果一致;但本研究中原始光谱特征波段位于近红外区域,有别于陈澜等[28-29]的研究结果,原因在于试验数据采集对象为长势差异明显的成熟期葡萄叶片,不同叶片样本色素含量差异小于细胞栅栏组织稳定性差异,故特征波段位于近红外区域而非可见光区域。光谱数据采集过程中会不可避免地引入环境噪声,微分变换是削弱噪声的有效方式之一[30],随着分数阶的增加,特征波段逐渐由近红外区域向红边靠拢,DSI的组合亦由近红外和绿光变为近红外和红边,这是因为当分数阶增加到一定程度时会导致不同SPAD值葡萄叶片微分光谱在可见光区域的一致降低,变相削弱该范围波段反射率对SPAD值的敏感性。本研究中特征波段、DSI与SPAD值的相关性随分数阶的递增呈先升后降趋势,二者相关系数分别在0.6和0.4阶达峰值,分别为0.777和0.7850可见分数阶微分变换虽弱化了可见光区域敏感性,但也提升了单因素模型建模变量反演SPAD值的潜力。

针对“维度灾难”,目前高光谱数据的降维手段主要分为特征提取和特征选择[31],前者通过特定算法重置响应变量的特征空间,对响应变量所有特征重投影形成不同于原有特征向量的新自变量,如主成分分析(PCA)和最大噪声比率(MNF)等,后者则基于特定条件从响应变量特征空间中选择部分特征作为自变量,如竞争自适应重加权采样(CARS)和无信息变量消除(UVE)等。本研究采用SPA算法对葡萄叶片原始及分数阶微分光谱降维,不同类型光谱被选择变量为5～16个,降维效果明显,但由于植被种类、采样设备及环境等的不同,SPA选择的波段位置及数量与刘燕德等[32]的研究结果不一致,后续的研究应着力于葡萄其他生育期和不同降维算法,归纳葡萄叶片SPAD值的关键反演波段与最佳降维方法。

本研究中多因素模型优于单因素模型,机器学习算法可对传统回归模型作进一步优化,与Wang等[33-35]的研究结果一致。机器学习算法本质上属于学习策略,通过模拟自然界中各种自然现象以解决实际问题,如神经网络算法通过模拟人类大脑神经元的信息反馈过程构建输入层、中间层和输出层以整体拟合目标特征,XGBoost则模拟树枝分叉的特性拆解目标的各项特征,将待监测对象转换为不同树枝特征之和,然而机器学习算法需要较大的时间成本确定各种参数和超参数,因此,学者们开发不同寻优算法配合以节约时间成本[36],未来应更多地寻优算法应用到XGBoost模型以明确最佳寻优算法。本研究在建立SPAD值反演模型时仅从数理统计的角度出发,并未考虑葡萄生长的生物学特征,所构建的反演模型是否适用于其他地区也还有待进一步测试;后续的研究应将数理统计模型与辐射传输模型融合,在其他区域开展相似试验,进一步提升反演模型鲁棒性与泛化性,同时模型自变量的选择应尽量贴近低空无人机遥感和卫星遥感,以期未来构建空-天-地一体化的葡萄叶片SPAD值高精度反演模型。

4 结 论

本文以贵州息烽葡萄种植基地成熟期葡萄作为研究对象,同步采集冠层叶片高光谱数据和SPAD值数据,基于不同分数阶微分光谱构建多个葡萄叶片SPAD值反演模型,结果表明:葡萄叶片SPAD值与原始光谱可见光区域反射率成反比,与近红外区域反射率成正比。随着分数阶的增加,单因素模型建模变量与SPAD值的相关性呈先升后降趋势,特征波段由近红外向红边靠近,DSI计算波段由近红外与绿光组合变更为近红外与红边组合,多因素模型建模变量逐渐剔除蓝光与近红外区域;除0.6阶和0.8阶外,其余分数阶微分光谱单因素模型建模变量均为DSI,其在特定分数阶下优于特征波段。多因素模型优于单因素模型,GA-XGBoost算法可进一步优化传统回归模型;所有构建的模型以SPA-GA-XGBoost模型精度最高,其建模与验证的R2分别为0.79和0.75,相应nRMSE分别为14.45%和15.54%。