基于机器学习的油青菜心水分胁迫研究

杨明欣，高 鹏，陈文彬，周 平，孙道宗,2，谢家兴,2，陆健强,2，王卫星,2

(1 华南农业大学 电子工程学院/人工智能学院，广东 广州 510642；2 广东省农情信息监测工程技术研究中心，广东 广州 510642)

油青菜心Brassicachinensisvar.parachinensis是一种重要的绿色蔬菜，广泛种植于广东省及全国多个城市，其叶片为油绿色，形状尖细，风味甜爽，菜苔油绿有光泽[1-2]。油青菜心对土壤水肥具有较高的要求，在生长过程中需要保证足够的水分供应，以满足其正常生长需求[3]。在持续缺水状态下油青菜心的多项生理活动受到影响，如叶绿素含量下降、叶片含水量下降等[4]。通过监测油青菜心的水分状况，采取高效的水分灌溉策略可保证油青菜心的品质，同时达到节约水资源、提高水资源利用效率的效果。

当作物缺水时冠层温度升高，基于冠气温差的水分胁迫指数能实时反映作物的水分亏缺状态[5-9]。水分胁迫指数 (Crop water stress index，CWSI)是反映作物水分状况的无量纲因子，范围在[0，1]之间，0表示作物无水分胁迫或者灌溉充分，1表示作物无蒸腾作用或严重胁迫。国内外研究在使用CWSI经验模型时，上下限方程确定方法并不统一[10]。Jones等[11]使用涂抹凡士林和水的干、湿参考面温度计算葡萄的CWSI；但是使用干、湿参考面需考虑位置对计算模型的影响[12]。王卫星等[13]采用水汽压差为0和6时对应的温度研究柳叶菜心的水分状况；但是水汽压差随着地域和季节变化，意味着基准线也会变化[14]。张立元等[15]、Agam等[16]使用空气温度+5 ℃作为上限研究玉米、橄榄树等作物的水分变化规律，Kumar等[17]、Khorsandi等[18]使用空气温度+2、+3 ℃作为上限研究芥菜、芝麻的CWSI，避免了建立上限方程；虽然能构建CWSI，但是不少研究直接套用经验值，没有结合作物本身及环境情况进行分析。屈振江等[19]研究苹果冠层温度的变化特征时发现，冠层温度峰值出现的时间、变化趋势与空气温度一致。针对相似的数据，具有噪声的基于密度的空间聚类(Density-based spatial clustering of applications with noise，DBSCAN)算法可将其划分到集合中，以簇中心代表集合的数据[20]。因此，本研究拟结合DBSCAN聚类算法，在给定的空气温度范围内探讨油青菜心冠气温差上限的固定值。

谢慧婷[21]、Matese等[22]研究生菜、葡萄藤的水分状况时发现，在缺水状态下光合作用速率呈下降趋势，而CWSI与光合作用速率呈负相关的关系。利用光合作用测量系统可直接获取作物的光合作用速率[23-24]，通过叶室夹住被测叶片形成固定被测空间并取样实现数据的自动采集，但叶室对叶片有一定程度的损坏。通过其他容易获取的数据模拟作物的光合作用速率，能有效解决损坏叶片的问题，刘煦[25]使用有效辐射、相对湿度等环境数据，基于机器学习方法预测了林下参的光合作用速率；陈硕博[26]、陈俊英等[27]利用棉花冠层的光谱数据，以线性回归、机器学习多种方法反演了光合参数。机器学习方法已被广泛应用于农业领域，宋飞扬等[28]、Botula等[29]使用最邻近节点(k-Nearest Neighbor，KNN)算法筛选特征、模拟土壤持水量，提高了模型的预测精度；Mohammadi等[30]使用支持向量回归 (Support vector regression，SVR)算法模拟每日参考蒸发量，结果表明预测值能较好地拟合真实值；Wang等[31]利用土壤回声结合极端梯度提升 (Extreme gradient boosting，XGBoost)算法预测土壤pH，达到预测精度高、误差低的效果；白婷等[32]应用光谱数据和随机森林(Random forest，RF)算法获得较高精度的土壤有机质估测值。本研究拟建立CWSI与光合作用速率的关系，以CWSI预测光合作用速率，探究不同胁迫程度下油青菜心光合作用速率的变化规律。在油青菜心从四叶一心到植株现蕾的营养生长阶段(Vegetative growth stage，V期)和从植株现蕾到菜苔高度与苔叶先端齐平的生殖生长阶段(Reproductive growth stage，R期)进行不同水分处理，计算CWSI，运用DBSCAN聚类分析方法选取作为冠气温差上限的固定值，采用 KNN、SVR、XGBoost、RF 这 4种机器学习方法预测光合作用速率。

1 材料与方法

1.1 试验材料

以油青菜心为试验材料，试验于2020年11—12月在华南农业大学工程学院进行。11月初将种子播种于花盆中，花盆高度15 cm，上口径22 cm，每盆装适量的土壤。待菜心出芽长出第1片真叶后，将相同长势的菜心移到相同规格的花盆中，每盆2株，间距为10 cm，待幼苗长到四叶一心时期开始进行不同水分处理，如图1所示。

图 1 油青菜心盆栽图Fig.1 The pot cultivation picture of Brassica chinensis

试验前每个花盆灌溉充足水分，静置1 h，每个花盆取表层0～10 cm的土壤，称质量(m)，放入烘箱在105 ℃条件下烘干6 h，待土壤冷却后再称质量(md)，按公式(1)计算该花盆土壤的田间持水量(w)[33]：

以平均值作为试验采用的土壤田间持水量，计算结果为32.2%。

1.2 水分处理及数据测定方法

根据本试验土壤的田间持水量(32.2%)，将试验对象油青菜心的水分胁迫梯度设置为5组，分别用T1～T5表示。T1的田间持水量为32.2%，以T1作为对照，T2～T5的田间持水量分别为85%T1、70%T1、55%T1 和 40%T1。

数据采集方法：选取菜心冠层顶部完全展开的、并能获取充足的阳光的叶片作为测量目标。数据采集时间为每天10:00—16:00，每30 min采集1次。使用手持式红外测温仪(型号Raytek ST18，雷泰公司产品)进行冠层温度测量，该设备的光学分辨率为12︰1，发射率为0.95，光谱响应范围8～14 μm，测温范围–20～500 ℃。测量时，保持红外测温仪与测量目标的距离在叶片大小的12倍以内。每个目标点测量3次，取平均值。使用手持式温湿度计(型号Aicevoos W8，艾沃斯公司产品)进行空气温度、相对湿度的测量，测量时保持设备与测量目标的距离为10 cm，1 min内测量3次后取平均值。使用土壤水分传感器(型号RS485，鹰都公司产品)进行土壤含水量测量，测量时将设备埋在2个目标叶片中间，与土壤表面的距离为10 cm。使用光合作用测定仪(型号SYS-GH30D，塞亚斯公司产品)进行光合作用速率、光合有效辐射的测量，该仪器基于快速准确的红外线CO2气体分析仪法，测定量程为 0～3 000 mg/kg，精度为 3 mg/kg，测量时使用叶室夹住叶片，保持1 min，测量5次后取平均值作为该目标点的实际光合参数。

受试验期间天气影响，在保证数据充足的前提下，本文剔除了阴雨天气(12月13—20日)测量的数据。

1.3 水分胁迫指数（CWSI）计算

CWSI的计算如以下公式所示[15]：

式中：θc为作物冠层温度；θa为空气温度；(θc−θa)ll为下限方程或无水分胁迫基准线，是作物在无水分胁迫时或充分灌溉下的冠气温差；(θc−θa)ul为上限方程或无蒸腾作用基准线，是作物在无蒸腾作用时、气孔关闭状态下的冠气温差；单位均为℃。

(θc−θa)ul和 (θc−θa)ll的计算如以下公式所示：

式中：A和B为回归系数；RH为空气相对湿度；VPD和VPG分别为饱和水汽压差和饱和水汽压差梯度。

根据不同的CWSI经验模型，分别使用公式(3)计算 (θc−θa)ul和将固定值作为 (θc−θa)ul，固定值由DBSC AN算法得到。

1.4 基于密度的聚类算法(DBSCAN)

DBSCAN算法能找到任意形状、集中分布的簇[34]。在给定的数据集中，DBSCAN算法将所有对象标记为未访问，随机选择未访问的对象p并标记为已访问，通过检查p的ε-(ε>0)领域内包含的MinPts个对象判断其属性，其中ε为半径参数，MinPts为领域密度阈值。若为核心点，将p的ε-领域中的所有对象添加到集合C中，直到C不再扩展，迭代停止。

距离指标常用欧式距离和曼哈顿距离，特征向量Xi、Xj之间的欧式距离、曼哈顿距离计算如公式(7)、(8)：

式中：dij表示第i、j个特征向量Xi、Xj的距离；k表示特征向量的维度；xi(m)、xj(m)分别为Xi、Xj在第m维的值，m的取值范围为1，2，…，k。

DBSCAN不需要预先设定簇数目，由算法自动决定，但需要设定半径参数ε和领域密度阈值MinPts。在本研究中，以欧式距离作为度量，ε为0.44，MinPts为 4。

1.5 最邻近节点算法(KNN)

KNN是数据挖掘算法中最简单的一种，精度高且对异常值不敏感[35]。对于要预测的点(xi,yi)，KNN在一系列样本坐标中选择k个离xi最近的样本坐标，对其y值求平均，结果为KNN模型的预测值。

预先设定的值包括k值和k个邻近点的权重，在本研究中，以曼哈顿距离作为度量，k为3，各邻近点的权重一致。

1.6 支持向量回归(SVR)

SVR被广泛应用于处理模式识别问题，泛化错误率低[36]。该算法结合了核函数和线性回归，让训练集中的每个点(xi,yi)尽量拟合到一个线性模型yi=ωx+b。将最优超平面记作 ωx+b=0，样本x到最优超平面的距离为r，通过引入松弛变量和惩罚系数求最小距离，将最优超平面问题转化为最优化问题：

最终支持向量回归模型如公式(13)所示：

1.7 极端梯度提升法(XGBoost)

XGBoost是Chen等[37]开发的一个开源机器学习项目，处理速度快，高度灵活，能更好地控制过拟合。该算法在梯度提升决策树基础上对损失函数进行二阶泰勒展开并添加正则项，通过不断形成新决策树拟合之前预测的残差，从而减少预测值与真实值残差。其目标函数L由损失函数l和正则项 Ω组成，如公式 (14)、(15)：

式中：L(∅)t表示第t次迭代的目标函数；yi为目标值；表示前t−1次迭代的预测值；l为目标值与预测值的平方差；ft(xi)为第t次迭代产生的新模型；Ω(fk)表示第t次迭代的模型的正则项；γ和λ表示正则项系数；T表示该模型的叶结点个数。

对公式(14)使用泰勒展开，可得：

式中：gi表示样本xi的一阶导数；hi表 示样本xi的二阶导数；ωj表示第j个叶子结点的输出值；Ij表示第j个叶子结点值的样本子集。对 ωj求导并且令导函数等于0，可求得使目标函数达到最小值的 ωj。在本研究中，学习步长为0.09，最大的树深度为4，弱学习器数量为70，正则化系数γ和λ为0.001和1。

1.8 随机森林算法(RF)

RF最早由Breiman[38]于2001年提出，该算法基于分类树，是Bagging算法的优化，具有运算速度快、稳定性好等优势，在处理大数据集时预测精度高。RF从总体数据集中采用多次、随机的方法取一部分样本构成样本簇，在每一个新生成的样本簇中训练回归决策树，组合每一棵决策树的输出结果加权平均。

对于给定的训练数据集，D={(x1,y1),(x2,y2),···,(xi,yi)}，其中xi为输入实例(特征向量)，yi为目标向量，i=1，2，…，n，n为样本容量。每次划分特征空间，逐一考察当前集合中所有特征的所有取值，根据平方误差最小化准则选择其中最优的一个作为切分点。对于训练集中第j个特征变量xj和取值s，作为切分变量和切分点，并定义2个区域，为找出最优j和s，对以下公式求解：

式中：c1、c2为划分后2个区域内固定的输出值；利用选定划分区域内yi求平均得到相应的输出值，直到生成完整的决策树。

1.9 模型评估

本研究运用SPSS软件，采用误差分析评估机器学习模型的预测效果，利用最小显著极差法分析水分处理间的显著水平。误差分析包括相关系数(Correlation coefficient，R2)、平均绝对误差 (Mean absolute error，MAE)、均方根误差 (Root mean square error，RMSE)，如公式 (20)～(22)所示：

式中：E为数学期望；f表示预测值；y表示实际值；i为第i个(i≤n)数据，n为数据量。

最小显著极差 (Least significant difference，LSD)法，对任何2个i、j处理平均数间的均数，若其绝对值≥ LSDα，则为在α的水平上差异显著，LSDα的计算如公式(23)所示：

式中：d fe为误差自由度；tα为误差自由度下的临界值；MSe为误差均方；n为各处理内的重复数。

2 结果与分析

2.1 充分灌溉下冠气温差与饱和水汽压差的关系分析

在无云晴天采集充分灌溉T1处理的冠层温度、空气温湿度数据，建立CWSI经验模型的冠气温差下限。图2a、2b分别为营养生长阶段(V期，11月27日—12月12日)生殖生长阶段(R期，12月21日—12月31日)的冠气温差与饱和水汽压差的散点图，通过回归方程拟合得到冠气温差的下限方程，方程如下所示：

图 2 冠气温差与饱和水汽压差(VPD)的关系Fig.2 Relationship between canopy and air temperature difference and vapor pressure deficit (VPD)

由图2可以看出，冠气温差与饱和水汽压差的拟合结果较好，两者显著相关(P<0.01)，均方根误差分别为1.54(V期)和1.07(R期)。在V期和R期冠气温差的下限有一定差异，R期的冠气温差下限比V期高，因此在建立油青菜心CWSI经验模型时应针对不同生长期建立相应的冠气温差下限。

根据所拟合的冠气温差下限方程和公式（3）计算CWSI，并绘制CWSI每日变化曲线，如图3所示。随着田间持水量降低，水分胁迫程度加深，CWSI明显增加。处理间的CWSI平均涨幅为0.13，具有显著差异(P<0.01)，表明CWSI能较好地反映油青菜心的水分状况。T1处理为充分灌溉，其每日平均 CWSI在 [0, 0.2]波动，T2、T3 和 T4 处理的每日平均 CWSI分别是 [0.2, 0.4]、[0.2, 0.5]和[0.3, 0.6]，受胁迫最多的T5处理的每日平均CWSI在[0.4, 0.8]波动，变化范围较大，表明水分胁迫程度越大，油青菜心越容易受到气候的影响。处理开始时变化曲线密集，随着胁迫时长增加，变化曲线无重叠或交点，表明水分的长期缺失对油青菜心的生理活动造成了一定影响，使得曲线波动范围缩小。

图 3 5种田间持水量处理的水分胁迫指数(CWSI)变化曲线Fig.3 Crop water stress index(CWSI) change curves under five field water holding capacities

2.2 不同生长期的冠气温差上限分析

根据公式(3)获取的冠气温差上限及对应的空气温度，运用DBSCAN算法对油青菜心进行聚类，结果如图4所示。由图4可以看出，随着空气温度上升，2个时期的油青菜心冠气温差上限均呈现上升的趋势，其中R期的冠气温差上限比V期高，主要分布在[4，5] ℃，V期的冠气温差主要分布在[3，4] ℃。由于试验期间空气温度变化较大，算法自动聚类的簇共有6个，且R期的空气温度变化比V期更大，因此R期有4个簇，V期有2个簇，表明空气温度对DBSCAN的聚类效果影响较大。跨度越大则簇越多，增加了选取固定值的难度。数据集主要分布在[20，30] ℃，随着数据的稠密程度上升，簇内包含的数据增多，其中2个大簇的质心分别为3.4 ℃(类别为 2，116)和 4.2 ℃(类别为 1，164)，分别对应V期和R期。此外，聚类结果显示有小部分噪声点(类别为−1，是离群点)，大多分布在空气温度区间边缘，剔除该部分数据并不影响聚类的效果。

图 4 基于DBSCAN的冠气温差上限聚类Fig.4 Clusters of the upper limit of canopy temperature minus air temperature by DBSCAN

为验证聚类的效果，根据无蒸腾作用基线获取的上限，分别采用空气温度+2.0、+5.0 ℃与聚类结果计算T1、T5处理的水分胁迫指数，其误差分析如表1所示。由表1可以看出，在2个阶段、2个处理中，空气温度+2.0 ℃的CWSI平均值偏高，空气温度+5.0 ℃的CWSI平均值偏低，在充分灌溉T1处理的条件下与无蒸腾作用基线的CWSI均无显著差异，当田间持水量为40%T1(T5处理)时与无蒸腾作用基线的CWSI均有显著差异(P<0.05)。空气温度+3.4、+4.2 ℃的CWSI均位于空气温度+2.0、+5.0 ℃之间，与无蒸腾作用基线的CWSI显著相关(R2=0.99)。不同生长期、同温度、同处理的CWSI存在显著差异(P<0.05)，当V期和R期采用同一个温度上限时，4.2 ℃将造成V期的CWSI偏小，3.4 ℃将造成R期的CWSI偏大，因此应结合作物的生长期选取不同的固定值。与充分灌溉T1处理相比，采用固定温度计算T5处理的CWSI与无蒸腾作用基线的平均值差值差异较大，V期2.0、3.4和5.0 ℃的涨幅分别为0.152、0.003和 0.104，R期 2.0、4.2和 5.0 ℃ 的涨幅分别为0.350、0.019和 0.044。拟合结果表明，3.4和 4.2 ℃可分别作为油青菜心在V期和R期的冠气温差上限，本研究将以2个温度计算的CWSI作为机器学习模型的输入向量之一。

表 1 不同冠气温差上限的水分胁迫指数(CWSI)误差分析Table 1 Error statistics of crop water stress index (CWSI) under different upper limits of canopy and air temperature difference

2.3 不同水分处理下的光合作用速率日变化分析

为研究油青菜心CWSI与光合作用速率的关系，绘制了光合作用速率的日变化曲线(图5)。由图5可以看出，光合作用速率呈现先增加而后减少的变化趋势，这是由于11:30时左右的太阳辐射较强、温度较高，有利于作物进行光合作用。在V期11:30时以后光合作用速率的下降趋势较明显，通过基于LSD法的多重比较分析，各个时刻的光合速率平均值差异均达0.05的显著水平，在R期光合作用速率的下降趋势较平缓。在V期、R期中，CWSI与光合作用速率是负相关的关系，充分灌溉处理T1的光合作用速率最高，当水分胁迫程度加深时，菜心叶面温度上升，光合作用速率下降，不同水分胁迫处理的光合作用速率在0.05水平具有显著差异(LSD法)。在V期，11:30时5组水分胁迫处理的光合作用速率差距最大，在15:00时差距最小。

图 5 不同水分胁迫处理的油青菜心光合作用速率日变化Fig.5 Diurnal variations of photosynthetic rate of Brassica chinensis in different water stress treatments

2.4 基于机器学习的光合作用速率预测分析

本试验运用 KNN、SVR、XGBoost、RF 这 4种模型，以CWSI和光合有效辐射作为输入向量预测油青菜心的光合作用速率。光合作用速率测量值与预测值的散点图、拟合方程如图6，数据经过标准化处理，预测模型的误差分析如表2。

图 6 4种模型的光合作用速率测量值与预测值的散点图Fig.6 The scatter plots of predicted and measured photosynthetic rates based on four models

表 2 4种预测模型的误差分析Table 2 Error statistics of four predicted models

由图6和表2可以看出，XGBoost模型的预测值与测量值的相关系数最高，其次是SVR、KNN和RF模型，4种预测模型的相关系数均大于0.85，拟合方程的截距及斜率均没有显著差异，表明这4种机器学习均可以预测油青菜心的光合作用速率。由于光合作用速率的整个数据集分布并不均匀，随机划分的测试集在各区间的分布与整个数据集基本一致，主要集中在[0.5,1.5]，当光合速率>2时预测值均偏大，表明原始数据的分布对这4种模型的预测效果有一定影响。集成学习方法XGBoost、RF在建模时往往需要规模大的数据集，本试验用于预测的数据集共有498组输入向量，XGBoost、RF的预测结果容易出现欠拟合的情况，无法体现提升树的优势。当数据集为原来的三分之一时，预测精度均下降，4种模型建模的相关系数分别为0.836(KNN)、0.868(SVR)、0.871(XGBoost)、0.837(RF)，下降幅度最大的是KNN和RF，表明SVR在小样本回归中比RF、KNN表现好。当数据集较小时，SVR和XGBoost均可通过调整参数提高预测效果，如SVR减少惩罚系数和松弛变量，XGBoost提高学习速率、降低最大生成树的数目，但KNN和RF并没有显著提高。XGBoost和RF都属于集成学习，但XGBoost使用了一阶导数和二阶导数，使得预测值与实测值的拟合残差减少，又借鉴了RF列抽样的做法，降低过拟合和减少计算，大大提升效率，因此预测效果最好。

3 结论

本研究通过5组水分处理的对比分析揭示了油青菜心在营养生长阶段和生殖生长阶段的水分胁迫变化规律，运用DBSCAN算法对油青菜心冠气温差进行聚类，探讨了使用固定值作为上限的可行性，采用KNN、SVR、XGBoost和RF算法预测了油青菜心的光合作用速率。结果表明，CWSI随田间持水量减少而升高，说明CWSI可以很好地识别油青菜心的水分胁迫；聚类结果显示，3.4和4.2 ℃与经验公式计算的CWSI具有显著的相关性，说明采用固定值计算油青菜心CWSI具有可行性，且2个阶段的温度存在差异，表明应结合作物的生长期进行聚类分析；4种机器学习方法均取得较好的预测效果，XGBoost模型的预测值与实测值拟合效果最好，其次是SVR、KNN和RF，说明机器学习方法预测光合作用速率具有可行性，且实现了光合作用速率的快速检测。本研究只建立了油青菜心在11—12月的水分胁迫指数模型，由于经验方程受天气、季节变化的影响，下一步试验将采集全年不同月份的数据，并考虑太阳辐射、风速的影响。在预测光合作用效率时，数据尚不够充分，作物种类单一，下一步研究将扩充数据集，并期望能在多个作物上应用。