基于高光谱和数据挖掘的油菜植株含水率定量监测模型

潘 月， 曹宏鑫， 齐家国， 吴 菲，3， 韩旭杰， 丁昊迪， 葛道阔， 张玲玲， 张伟欣， 张文宇

(1.南京农业大学农学院/亚洲农业研究中心，江苏南京210095；2.江苏省农业科学院农业信息研究所，江苏南京210014；3.扬州大学农学院，江苏扬州225009)

油菜是中国种植面积最大的油料作物之一和国产植物油第一大油源，2020年中国油菜种植面积约6.80×106hm2，油菜籽产量1.400×107t，仅次于加拿大，居世界第二，发展油菜生产是保障中国食用油供给安全的根本途径[1]。油菜是需水量较大的作物[2-4]，植株水分状况不仅影响油菜生长[5]，还会影响其产量、氮素吸收能力以及油菜籽的品质[6-9]。而伴随水资源日益紧张，传统的“丰水高产”农业生产模式逐渐向现代“节水高产”模式转变，水分动态定量诊断与调控是实现这个转变的关键技术，其基础则是植株水分状况的定量监测。高光谱手段能够采集精细的光谱信息，其在监测植物水分状况中的应用十分广泛。运用遥感技术监测植物中的水分，学者更倾向于选择冠层和叶片指示水分状况[10]。随着传感传输技术的快速发展，遥感技术能够改变以往费时费力、破坏性大的诊断方法，实现植物水分实时监测、精确诊断。在单叶尺度方面，杨玉清等[11]用偏最小二乘法(Partial least squares, PLS)及逐步多元线性回归(Stepwise multiple linear regression, SMLR)法对玉米叶片水分含量进行建模分析，验证集决定系数分别为0.975和0.980；Krishna等[12]研究了10个水稻基因型在不同水分胁迫下的表现，通过光谱指数、多元技术和神经网络筛选最优波段，建立预测模型，结果表明,偏最小二乘回归-多元线性回归(PLSR-MLR)是预测叶片含水率(RWC)的最优模型，验证R2为0.97。在冠层尺度方面，刘晓静等[13]的研究结果表明，在拔节期、抽穗期、开花期、灌浆前期和灌浆后期，冬小麦叶片相对含水量与归一化差值植被指数(NDVI)、优化土壤调整植被指数(OSAVI)、比值/归一化植被指数(R/ND)、冠气温差(TDC-a)和冠气温差(TDC-a)拟合效果较好，决定系数分别为0.842、0.884、0.831、0.864和0.945。精确监测作物水分状况的主要步骤包括分析光谱响应规律、筛选模型变量、构建监测模型和检验监测模型。在油菜水分状况监测方面，张晓东等[14-15]基于高光谱数据，采用区间分段逐步回归法定量分析油菜含水率，模型预测值与实测值间的相关系数为0.87，后续与多光谱图像和冠层温度多信息融合，将相关系数提高至0.93；仝春艳等[16]提出2种改进型角度指数ANI1 450和ASI1 450，在不区分苗期、蕾薹期的情况下估算叶片等效水厚度效果最好，R2均达到0.832；张君等[17]的研究结果表明，正交信号校正(OSC)+主成分分析(PCA)+SVR模型估测叶片含水率效果最好，训练集决定系数(Rc2)与测试集决定系数(Rv2)分别高达0.901和0.857。可见，已有研究多能够在单一生育期实现精确监测，而对于整个生育期的普适性有待提高，目前光谱分析方法在作物水分监测方面运用得较少，高光谱信息有待深度挖掘，机器学习方法在模型构建上应用得较少，且局限于实验室条件。本研究旨在基于不同品种、氮素与水分处理下的油菜池栽田间试验，利用逐步回归分析法(Stepwise regression, SR)、连续投影算法(Successive projection algorithm, SPA)和竞争自适应加权算法(Competitive adaptive reweighted sampling, CARS)筛选最优波段组合，并利用减量精细采样法(Reduced precise sampling method, RPSM)筛选最优归一化差值光谱指数(NDSI)和比值光谱指数(RSI)，为监测模型提供变量，利用线性回归(Linear regression, LR)、BP神经网络(Back-propagation neural network, BPNN)和支持向量机回归(Support vector regression, SVR)方法构建油菜植株含水率(PWC)监测模型，以进一步提高油菜PWC监测精度，为油菜水分诊断调控提供依据。

1 材料与方法

1.1 试验设计

试验于2019-2020年和2020-2021年在江苏省农业科学院实验农场(32.03°N, 118.87°E)进行，前茬为休闲，两季试验的基本情况如表1所示。供试品种3个，分别为浙杂903(C1)、宁油22(C2)和宁杂1818(C3)，施肥水平分别为不施肥(N0)、施纯氮180 kg/hm2(N2)，移栽密度统一为1 hm21.2×105株(行距40 cm,株距20 cm)。土壤基础养分：有机质含量31.4 g/kg，全氮含量2.03 g/kg，速效磷含量20.3 mg/kg，速效钾含量139 mg/kg，pH 7.31。

2019-2020年供试品种为浙杂903(C1)和宁杂1818(C3)，2020-2021年供试品种为宁油22(C2)和宁杂1818(C3)。采用裂区设计，主区为施肥，设为不施肥(N0)和施肥(N2，于移栽前施入有机肥60 kg/hm2，于抽薹期叶面喷硼15 kg/hm2，氮肥180 kg/hm2按基肥∶腊肥∶薹肥=5∶3∶2分配)处理，副区为水分与品种处理，水分处理设浇水2 min(约小区水层4.87 mm)(W1)、3 min(约小区水层7.31 mm)(W2)、4 min(约小区水层9.75 mm)(W3) 3个水平。2020-2021年浇水后不同处理0～10 cm耕层土壤体积含水率变化情况如图1所示，共12个处理，重复3次，随机排列。小区面积3.77 m×2.77 m≈10.44 m2，行距0.4 m, 株距20.0 cm。分别于2019年10月1日、2020年10月7日播种，分别于2019年11月2日、2020年11月11日移栽，其他栽培措施同高产大田管理。

表1 2个生长季试验的基本情况

图1 2020-2021年不同处理0～10 cm耕层土壤体积含水率变化Fig.1 Variation of volume water content of 0-10 cm topsoil under different treatments in 2020-2021

1.2 测定项目与方法

1.2.1 冠层光谱反射率 选择晴朗无云无风天气，分别于各个生育时期当天10:00-14:00在每个小区选3个代表植株并挂牌标记序号，每个植株测定5个数据，取平均值代表该小区的光谱数据。用美国ASD(Analytical Spectral Device)公司的ASD Handheld 2型手持式光谱分析仪距油菜冠层上方50 cm测定冠层光谱反射率，光谱仪视角为25°，测定波段为325～1 075 nm，光谱采样间隔为1.0 nm，光谱分辨率为3 nm，每次测量前用白板校正。

1.2.2 油菜植株含水率(PWC) 基本与光谱测定同步，分别在2019-2020年越冬期(移栽后77 d)、2020-2021年苗期(移栽后7 d、19 d、34 d、53 d)、越冬期(移栽后69 d)、终花期(移栽后149 d)取样，每小区取上述挂牌标记植株3株。鲜质量利用称质量法，样本称量鲜质量后，恒温105 ℃杀青30 min，之后恒温80 ℃烘干至恒质量，分别称量样本干质量。计算植株含水率，PWC=(PWF-PWD)/PWF×100%，式中PWF为植株地上部鲜质量总和，PWD为植株地上部干质量总和。取3株平均值代表该小区油菜植株含水率。

1.3 数据处理与分析

1.3.1 监测模型变量的筛选 本研究主要选用特征波段反射率和光谱指数为模型变量，运用Matlab软件计算，用Origin作图。

1.3.1.1 SR分析 SR分析以线性回归方程为基础，采用双检验原则，在建立“最优”回归方程的过程中按照贡献度由大到小的顺序，逐步引入对油菜PWC影响显著的波段，剔除影响不显著的波段，即每一步都进行F检验，保留F≥4时引入的波段，剔除F≤3.8时引入的波段，以筛选出最优波段组合。

1.3.1.2 SPA法 首先设定波长数(N)，然后根据不同波长数下的均方根误差(RMSE值)确定最佳的建模波长数量。每次向前循环选择时，都计算新引入的波长在剩余波长变量上的投影，保证它在选择的前一个波长变量正交子空间上具有最大投影值，从而最大限度地消除共线性信息[18]。

1.3.1.3 CARS法 首先利用蒙特卡罗采样法(Monte Carlo sampling, MSC)采样50次，运用指数衰减函数(Exponentially decreasing function, EDF)去掉权重较低的变量，运用自适应重加权采样技术(Adaptive reweighted sampling, ARS)筛选出偏最小二乘回归(Partial least regression, PLSR)模型回归系数绝对值大的变量，每次产生的新变量组合会再构建PLS回归模型并计算交叉验证均方差(RMSECV)，RMSECV最小的变量子集即为最优子集[19]。

1.3.1.4 RPSM法 采取RPSM法[20]，在680～1 075 nm波段范围内，首先每10 nm取光谱反射率均值，构建任意两波段组合的归一化光谱指数(NDSI)和比值光谱指数(RSI)，结合油菜植株含水率建立线性模型，计算决定系数，绘制等高线图以确定核心波段。对核心波段精细采样，以1 nm为间隔构建光谱指数。重复上述操作，选择决定系数最高的光谱指数作为模型变量。

NDSI=(Rλ1-Rλ2)/(Rλ1+Rλ2)

(1)

RSI=Rλ1/Rλ2

(2)

其中，Rλ1、Rλ2为680～1 075 nm范围内任意波段的光谱反射率。

1.3.2 监测模型的构建

1.3.2.1 线性回归 以筛选的光谱指数或敏感波段反射率作为自变量，油菜PWC作为因变量，运用最小二乘逼近拟合它们的定量关系，构建一元或多元线性回归方程。

1.3.2.2 BP神经网络 本研究定义隐藏层神经元数量为10，使用Levenberg-Marquardt反向传播算法进行训练。它是使用最广泛的非线性最小二乘算法，运用模型函数对估测数据在其领域内作线性近似分析，利用泰勒公式，忽略二阶以上导数项，能够使BP神经网络具有良好的非线性映射能力、容错性和较强的泛化能力。

1.3.2.3 支持向量机回归(SVR) 本研究使用支持向量机算法库LIBSVM[21]，SVM模型类型为epsilon-SVR，选用RBF核函数，损失函数P值为0.01，惩罚参数(c)和核函数参数(g)由穷举法确定。

1.3.3 监测模型的检验 采用常用的决定系数(R2)、均方根误差(RMSE)、平均绝对误差(MAE)和平均绝对误差占观测值的比值(Dap)对模型的精度进行评估，R2为监测值与实际值的决定系数，RMSE为监测的误差变幅，R2越高，RMSE、MAE和Dap越低，表明模型监测效果越理想。计算公式如下：

(3)

(4)

(5)

其中,n为模型测试的样本数，Pi为模型的估测值，Qi为试验的实测值。

2 结果与分析

2.1 油菜植株含水率与冠层光谱反射率的相关性

2个生长季油菜PWC与冠层光谱反射率的相关系数随波长的变化如图2所示。油菜PWC与冠层光谱反射率在325～1 075 nm波段大部分呈现极显著负相关性，即随着PWC增大冠层反射率降低，符合水分吸收太阳辐射使得光谱反射率降低的遥感原理，可见油菜冠层光谱反射率对PWC具有良好的指示作用。由于红边区域的680～780 nm反射率对水分有间接响应，近红外波段反射率对水分有直接响应，故选取红边波段和近红外波段进行后续分析。

图2 油菜植株含水率与冠层光谱反射率的相关系数随波长的变化Fig.2 Variation of the correlation coefficient between rapeseed plant water content and canopy spectral reflectance as the wavelength changes

2.2 油菜植株含水率监测模型变量的选择

2.2.1 基于SR分析的特征波段反射率 针对油菜PWC利用全波段数据进行SR分析，结果如表2所示。共建立了3个回归模型，模型R2随着波长数的增多逐渐增大，标准误差逐渐减小，所有模型均达到0.01显著性水平。第3个模型引入1 071 nm、986 nm和730 nm波长的光谱反射率，构建的回归模型R2为0.824，拟合度最佳，因此，将模型3作为最优波段组合用于构建油菜PWC监测模型。

表2 基于逐步回归分析筛选的油菜植株含水率(PWC)敏感波段

2.2.2 基于SPA的特征波段反射率 对油菜PWC运行SPA算法，其RMSE值随模型变量数(波长数)的变化趋势如图3所示，先上升后下降，波长数量大于7个时，RMSE值的变化不再显著，此时RMSE为1.33%。因此选取7个波长作为特征波长组合，分别为686 nm、695 nm、707 nm、746 nm、964 nm、1 065 nm和1 069 nm，其中4个位于可见光波段，3个位于近红外波段。

RMSE：均方根误差。图3 基于SPA算法的油菜植株含水率(PWC)的均方根误差(RMSE)随模型变量数变化趋势Fig.3 Variation trend of root mean square error (RMSE) for rapeseed plant water content (PWC) with number of model variables based on SPA algorithm

2.2.3 基于CARS的特征波段反射率 对油菜PWC运行CARS算法的筛选过程如图4所示。随着采样次数增加，在粗选过程中RMSE值迅速下降，然后下降速度逐渐变缓,转换到精选过程。剔除了与油菜PWC无关的波长，在35次采样后RMSECV逐渐增大，说明在筛选过程中丢失了有关信息。因此选择第35次采样筛选的10个波长为最优波长组合，分别为694 nm、695 nm、696 nm、863 nm、864 nm、893 nm、973 nm、986 nm、1 050 nm和1 071 nm，其中，3个波长位于可见光波段且十分集中，7个波长位于近红外波段。

2.2.4 基于RPSM的光谱指数比较 基于油菜PWC与任意两波段组合的NDSI决定系数如图5所示。初步采样全波段循环的决定系数如图5a所示，NDSI决定系数最大值均为0.716，近红外波段880～1 000 nm的部分决定系数超过0.500，因此选取880～1 000 nm的光谱反射率进行精细采样，结果如图5b所示，NDSI(R981,R894)为监测油菜PWC所得最优NDSI，决定系数为0.787。RSI决定系数情况与NDSI相似，RSI(R981,R894)为最优RSI，决定系数为0.786，略低于NDSI(R981,R894)，选其作为变量构建油菜PWC监测模型。

a：模型交叉验证均方根误差变化趋势；b：各波长回归系数随采样次数变化趋势。RMSECV:交叉验证均方差。图4 基于CARS筛选植株含水率(PWC)监测模型变量的过程Fig.4 Process of selecting variables of plant water content (PWC) monitoring model based on CARS

a：在680～1 075 nm范围内每隔10 nm采样；b：在880～1 000 nm范围内每隔1 nm采样。图5 油菜植株含水率(PWC)与任意两波段组合的归一化差值光谱指数(NDSI)的线性模型决定系数等高线Fig.5 Contour map of linear model determination coefficient (R2) of rapeseed plant water content (PWC) and normalized difference spectral index (NDSI) by any two band combinations

2.3 油菜植株含水率监测模型的构建及检验

2.3.1 线性回归模型 对于采集的60个油菜样本，42个用于线性回归模型构建，18个用于检验，模型效果如表3所示,建模R2由高到低依次为SPA、CARS、SR、RPSM-NDSI、RPSM-RSI。基于NDSI、RSI构建的模型建模和检验效果均相似，有较好的监测能力，但相比于SR、SPA、CARS筛选的敏感波段反射率略有不足，因为参与的波段数较少，丢失了一些有效信息。SPA-LR的建模R2和检验R2最高，分别为0.878和0.693，检验RMSE和Dap分别为1.623%和7.530%。

表3 油菜植株含水率(PWC)的线性回归模型及检验

建模R2为模型拟合度的决定系数，r为相关系数，R2为决定系数，RMSE为均方根误差，MAE为平均绝对误差，Dap为平均绝对误差占观测值的比值。**表示在0.01水平影响显著，r(0.01,16)=0.59 。SR:逐步回归分析；SPA：连续投影算法；CARS：竞争自适应加权算法；RPSM-NDSI：减量精细采样法-归一化差值光谱指数；RPSM-RSI：减量精细采样法-比值光谱指数。

2.3.2 BP神经网络模型 对于采集的60个油菜样本，42个用于训练BP神经网络模型，9个用于验证，9个用于测试，模型效果如表4所示，训练R2由高到低依次为SPA、RPSM-NDSI、RPSM-RSI、SR及CARS。以SPA-BP模型的训练R2最高，为0.964，SR-BP模型的验证R2最高，为0.977，SPA-BP模型的测试R2最高，为0.940，RMSE和Dap分别为1.836%和17.566%。

表4 油菜植株含水率(PWC)的BP神经网络模型及检验

2.3.3 SVR模型 对于采集的60个油菜样本，42个用于SVR模型训练，18个用于检验，模型效果如表5所示，训练R2由高到低依次为SPA、CARS、SR、RPSM-NDSI和RPSM-RSI。SPA-SVR模型的训练R2最高，为0.910，检验R2亦最高，为0.841，RMSE和Dap分别为1.227%和5.368%。

表5 油菜植株含水率(PWC)的SVR模型及检验

3 讨论

作物水分状况的光谱监测研究始于20世纪70年代，Thomas等[22]的研究结果表明，植物叶片的光谱反射率随叶片含水量下降而增加，1 450 nm和1 930 nm波段的光谱反射率和叶片含水量呈显著相关(P=0.01)。后续Holben等[23]、Jackson等[24]以及很多学者在不同植物上的研究结果都表明，主要水分吸收波段为0.97 μm、1.19 μm、1.45 μm、1.97 μm和2.70 μm[25]。高光谱能够更加精细地体现光谱的变化，因此除了主要水分吸收波段，亦存在其他能够反映作物水分差异的敏感波段，例如680～780 nm区域的红边光谱，PWC越低近红外波段反射率越高，红边区域的光谱反射曲线会更陡峭[26]，将这些波段加以利用，可以实现作物水分状况的精确监测，故本研究选取具备理论支撑的680～1 075 nm波段进行分析。然而高光谱数据量大，需要采取合适的手段提取关键波段，随着光谱分析技术的发展，全波段光谱快速分析成为可能。SR分析通过逐一引入光谱变量，剔除对监测目标影响不显著的变量，保留影响显著的变量，来筛选最优光谱变量集。SPA法是一种前向变量选择算法，能够使矢量空间共线性最小化，消除原始光谱数据中的冗余信息，该方法在叶绿素、生物量的高光谱定量监测中广泛应用[27-30]。CARS法基于“适者生存”的原则，以迭代和竞争的方式筛选最优波段组合。这些方法均能够在降低数据维度、筛选最优波段组合的同时充分挖掘数据、提取有效信息，从而监测并诊断作物水分状况[31-33]。光谱指数简单且具有较好的推广性，其中NDSI和RSI广泛用于获取植被生物物理特征属性[34]，减量精细采样法能够有效融合降采样法和精细采样法筛选最优光谱指数[20]。已有研究大多是运用一个方法筛选敏感波段，这些方法遵循的筛选原则不同，适配的建模方法不同，有必要和不同的建模方法组合，寻找筛选变量-建模效果最好的组合。

前述3个方法筛选的波段反射率变量更多、更分散，因此蕴含的信息更全面，在建模效果中即有体现，光谱指数构建模型的建模效果、检验效果低于全波段光谱分析。本研究旨在提高油菜水分状况的监测精度以及普适性，所用SR、SPA和CARS 3种全波段分析方法筛选的变量数量有所不同，但拟合精度最高，研究更侧重于比较不同变量筛选方法搭配不同建模方法构建的油菜PWC定量监测模型，今后的研究中可以采用更多的光谱分析方法，采集全波段有效信息，实现光谱数据的充分利用。

已有油菜水分状况的监测研究大多是运用线性或者非线性普通回归模型，在监测单一生育期作物水分状况时效果较好，而对全生育期的普适性有待提高。本研究引入BP神经网络、SVR进行非线性拟合，与线性回归模型相比较，模型拟合优度明显提高，能够为油菜水分诊断提供依据。本研究设置的水分梯度处理为常规油菜生存环境所面临的干旱、正常和渍水情况，对于特殊气候地区，有必要进一步扩大油菜PWC范围，以扩大模型的监测适应性。

4 结论

本研究基于不同品种、氮素与水分处理的油菜池栽田间试验，运用RPSM法、SR分析法、SPA法和CARS算法筛选监测模型变量，充分挖掘高光谱数据，构建线性回归模型、BP神经网络模型、SVR模型监测油菜PWC，模型的估测值和实测值均呈现极显著相关关系。其中SPA-LR模型估测值与实测值的检验R2为0.693，模型的RMSE为1.623%；SPA-BP模型估测值与实测值的检验R2为0.940，模型的RMSE为1.836%；SPA-SVR模型估测值与实测值的检验R2为0.841，模型的RMSE为1.227%。上述模型均能够实现油菜PWC的精确监测，以SPA-SVR模型的监测效果为最好，这可为油菜水分诊断提供依据。

致谢: 本研究所用试验材料宁油22和宁杂1818由江苏省农业科学院经济作物研究所张洁夫老师和浦惠明老师提供, 浙杂903由江苏省兴化市农业技术推广中心王洁站长提供,特致谢忱！