基于机器学习方法的杨树抗涝性评价和预测方法

技术领域

本发明属于机器学习和生物统计技术领域，具体涉及一种基于机器学习方法的杨树抗涝性评价和预测方法。

背景技术

洪涝胁迫不仅是指突发性或者永久性的土壤内涝，而是指所有不同程度下的土壤饱和，即杨树生长时地下水位远远高于杨树根部的土壤孔隙。这种条件下，就会导致土壤中的氧气水平下降，而乙烯碳和二氧化碳的浓度增加。杨树叶片的叶色素浓度发生改变，杨树不能正常进行光合作用。涝渍显著降低了杨树幼苗的苗高、地径和减弱光合作用能力等。涝渍灾害可进一步细分为涝害和渍害，二者对于植物生长的影响也不相同。涝害对植物幼苗的影响大于渍害，这是由于幼苗阶段整体抗涝机制还不成熟。洪涝胁迫条件下，植物的净光合速率、蒸腾速率和叶片气孔导度都会有所降低，并且随着土壤内涝时间延长而逐渐降低。洪涝灾害发生在土壤被水饱和时，导致植物根区处于厌氧状态，有氧呼吸受阻。全球气候变化使涝渍事件发生频率增加，造成相当大的农作物和森林植物损失。植物通过不定根生长、通气组织形成、能量代谢和植物激素信号传导等途径响应淹水胁迫。

内涝胁迫是杨树主要的非生物胁迫因素，通过改变杨树的生理代谢、延迟生长和发育的正常时间以及增加对疾病的易感性和感染程度来损害其正常生长。洪涝胁迫往往导致杨树只能在低氧甚至极端缺氧条件下生长，可直接诱发杨树内部大范围的伤害，如抑制有氧呼吸的正常进行、减少活性氧的产生和降低杨树光合速率。一般来说，杨树对这种低氧胁迫的反应生存能力，取决于它们对内涝胁迫的适应性和应对机制，而同一杨树不同品种个体之间，甚至不同基因型会表现出不同的耐涝性。

全球气候变化和洪涝灾害事件的变化息息相关，洪水通常会导致植物全部淹没或部分淹没，破坏当地土壤和地质条件，这对农业生产和生态系统都会重大损害。植物从种子萌发到营养和生殖生长都会受其影响，洪水胁迫导致植物在低氧条件下生长，从而对植物生命周期的发育期产生负面影响。在这种胁迫条件下，不同植物物种的生存策略是不同的，但是与洪水胁迫相关的感知途径在分子水平上是相似的。植物通过能量代谢、光合作用、呼吸和内源性植物激素生物合成信号传导的变化来响应洪水胁迫，由于有氧呼吸在洪水胁迫下受到抑制，能量代谢的减少进一步限制了植物的发育。

植物的生长通常会受到洪涝灾害的影响，而杨树具有生长周期短、耐涝性强的突出优点，是洪涝泛滥地区进行植树造林、改善生态系统重要树种。然而植物造林并不能完全抑制洪涝灾害的发生，同时洪涝灾害仍然会对这些地区的杨树生长和生存产生较大负面影响。因此，了解杨树如何感知和适应内涝或洪水覆盖的土壤深度大小的重要问题，以及赋予杨树耐涝性的相关机制具有重要意义，这将有助于筛选和培育具有更高耐涝性的杨树无性系品种。

目前对杨树的抗涝能力的预测研究尚不完善，尤其在对杨树抗涝能力相关地品种和特征的筛选实验方面，存在着实验周期过长、所需杨树品种和特征资料不足、筛选过程一般过于繁琐等问题。

发明内容

本发明的目的就是为了解决上述背景技术存在的不足，提供一种基于机器学习方法的杨树抗涝性评价和预测方法。

本发明采用的技术方案是：一种基于机器学习方法的杨树抗涝性评价和预测方法，包括以下步骤：

确定杨树抗涝性特征和抗涝性指标；

基于杨树抗涝性特征和抗涝性指标采用多种机器学习方法分别建立多个回归模型；

对多个回归模型分别进行模型评价，得出每个回归模型的GPI指数，根据多个GPI指数筛选出模型拟合效果较高的多个回归模型；

基于筛选的多个回归模型对每一个杨树抗涝性特征的重要性进行分析，得出杨树抗涝性特征的特征重要性排名，通过特征重要性排名评价和预测杨树的抗涝性。

进一步地，根据试验获取的若干初始杨树抗涝性特征，先基于专家进行初步筛选出若干个次要特征；然后通过递归特征消除法对若干个次要特征做进一步筛选，得到最终参与建模分析的若干杨树抗涝性特征。

进一步地，利用抗涝系数将杨树各生长指标进行加权组合形成抗涝性指标。

进一步地，所述抗涝性指标通过以下公式确定：

Pscore＝ω

其中，Pscore为抗涝性指标，ω

进一步地，还包括利用四分位法删除抗涝性指标的部分离群值。

进一步地，所述离群值的筛选公式为：

Q

其中，Pscore为抗涝性指标，Q

进一步地，所述多种机器学习方法包括XGboost、lightGBM、GradientBoost、DecisionTree、RandomForest和Adaboost 6种机器学习方法。

进一步地，所述GPI指数由评估模型性能的四个评价指标归一化得到。

进一步地，所述GPI指数通过以下公式确定：

其中，GPI

更进一步地，基于筛选的多个模型采用特征重要性和平均SHAP值两种方法对每一个杨树抗涝性特征的重要性进行分析。

本发明的有益效果是：

本发明利用苗高、地径和生物量的抗涝性系数构造抗涝性指标Pscore，将生物学问题转化为数学模型问题，不仅节省实验成本，而且能够包含更丰富的信息，增加了研究结果的代表性。

本发明通过机器学习方法的引入，使得杨树抗涝能力方面的分析有了更方便的手段和技术，对于接下来重要特征分析和杨树品种筛选都有很大改进。

本发明引入GPI指数对机器学习算法构建的预测回归模型进行排序，又利用特征重要性和平均SHAP值两种方法进行特征排序，筛选出了影响杨树抗涝性的重要特征，且利用这些重要特征进行预测，能有效杨树抗涝性评价和预测的精确度。

本发明提出的抗涝性指标以及相关回归模型在预测方面精确度和效率更高，能够分析出影响杨树抗涝性能力的重要特征及其特征排序，对进一步预测杨树抗涝性和选择强抗涝性的品种具有重要意义，也对开展抗涝性杨树品种种植提供了一定的参考价值。

本发明采用线性相关性热图和递归特征消除法(RFE)对杨树抗涝性相关的特征进行特征筛选，分析出和抗涝能力相关的重要性比较靠前的特征，不仅可以节约实验成本和实际种植成本，而且为接下来的杨树品种筛选和洪涝多发地区植树造林提供有益的指导。

附图说明

图1为本发明的流程图。

图2为本发明抗涝性指标Pscore与单个指标对比图。

图3为本发明基于Mean Shap Values对LightGBM模型特征排序图。

图4为本发明基于Permutation Importance对LightGBM模型特征排序图。

图5为本发明基于Permutation Importance对GradientBoost模型特征排序图。

图6为本发明基于Mean Shap_Values对GradientBoost模型特征排序图。

具体实施方式

下面结合附图对本发明的具体实施方式作进一步说明。在此需要说明的是，对于这些实施方式的说明用于帮助理解本发明，但并不构成对本发明的限定。此外，下面所描述的本发明各个实施方式中所涉及到的技术特征只要彼此之间未构成冲突就可以互相结合。

如图1所示，本发明提供一种基于机器学习方法的杨树抗涝性评价和预测方法，包括以下步骤：

确定杨树品种特征和抗涝性指标；

基于杨树品种特征和抗涝性指标采用多种机器学习方法分别建立多个回归模型；

对多个回归模型分别进行评价得出每个回归模型的GPI指数，根据多个GPI指数筛选出模型拟合效果高的多个模型；

基于筛选的多个模型对每一个杨树品种特征的重要性进行分析，得出杨树品种特征的特征重要性排名，通过特征重要性排名评价和预测杨树的抗涝性。

实施例

1、数据获取

为了更方便了解和研究杨树在淹水胁迫下的生长情况，本实施例利用淹水胁迫试验模拟杨树真实的涝渍生长环境。淹水试验设计中共有杨树品种21个，小叶杨在试验过程中死亡故其数据全部舍弃，不再参与后续的抗涝性研究。在本试验中，20个杨树品种分别设置4个实验组和对应的4个对照组，共有160种杨树材料。试验中将20个杨树无性系的1岁分枝切成约15厘米的插枝，有3或4个芽。在水中将准备好的插枝浸泡24小时后，再将插枝种植在幼苗花盆混合土壤中。试验中每天都能观察到杨树幼苗的具体生物形态变化，包括褪绿、叶片发黄和叶片脱落以及幼苗高度等。

淹水试验测定了杨树幼苗在第0天、第15天、第30天、第45天和第60天的幼苗高度、地径、生物量、光合作用和叶绿素荧光参数。叶片气体交换测量参数主要包括净光合速率(Pn)、细胞间二氧化碳分压(Ci_Pa)、细胞间二氧化碳浓度与环境二氧化碳浓度比值(Ci/Ca)和植物蒸腾速率(Tr)；以及周围环境条件参数(同一时间段中)，如叶片表面温度、光合活性辐射、相对大气温度以及相对空气湿度。在本试验获取数据过程中，最终测定了关于20个杨树品种幼苗在淹水胁迫试验前后的48个相关试验参数，包括光合作用相关参数、叶绿素荧光特征参数和环境变量参数。

2、数据预处理

实验获取的数据因为记录不完整或丢失而有部分缺失值，由于实验组和对照组都采取了四组实验，故缺失值一般采用每组的剩余值的均值进行填充处理。由于淹水导致死亡或者其他因素产生了部分异常值。根据实验组和对照组的其他组数据，对于每组的异常值选择直接删除的方法。为了检验杨树品种的抗涝能力，接下来对每个杨树品种的实验组数据取均值，并采用四分位法对极端值进行处理以保证后续的建模效果。

为了筛选抗涝指标，我们需要对生物量、地径和苗高的实验组和对照组数据进行显著性差异分析。实验获取的数据进行缺失值和异常值处理之后，首先利用正态检验Q-Q图进行分析：如果通过正态性检验，则对该指标实验组和对照组数据进行正态检验分析是否存在显著差异；对于不通过正态性检验，采用Mann-Withney-Wilcoxon检验分析指标是否存在显著差异。

3、确定杨树抗涝指标Pscore

获取上述实验数据后，本实施例利用抗涝系数将苗高、地径和生物量归一化，得到新的抗涝性指标Pscore，旨在充分衡量杨树在洪涝胁迫条件下生存情况。由于15天、30天、45天的杨树幼苗生长情况不太明显，且部分数据无法在每一阶段都进行充分测量和记录，故本实施例选择各杨树品种的苗高、地径和生物量的60天试验前后实验组和对照组样本数据变化情况作为衡量杨树抗涝性指标性能的重要参考。为了筛选出抗涝性指标，首先将苗高、地径和生物量标准化，标准化之后可以消除量纲和取值范围差距的影响，标准化公式分别如下所示：

其中，m、d、s分别代表苗高、地径和生物量在实验组的第0天和第60天变化量差异数值，mean(m)、mean(d)、mean(s)分别为实验组中苗高的均值、地径的均值和生物量的均值；std(m)、std(d)、std(s)分别为实验组中苗高的标准差、地径的标准差和生物量的标准差；Mscore、Dscore和Sscore分别代表苗高、地径和生物量标准化处理后的数值，标准化后三个指标取值范围是(-1,1)，处理之后更有利对三个指标进行定量分析研究。

为了得到新指标Pscore，可以利用抗涝系数将三个标准分数加权起来，抗涝系数是每个指标(苗高、地径和生物量)实验组样本总量与对照组样本总量的比值，抗涝系数的计算公式如下所示：

其中，m

利用抗涝系数将苗高、地径和生物量的标准化处理后的数值Mscore、Dscore和Sscore加权为一个新的抗涝性指标Pscore，其公式为：

Pscore＝ω

其中ω

在杨树指标单因素方差分析ANOVA(表1)中，将苗高、地径和生物量当作同一因素不同水平进行检验，检验结果显示F值是104，P值<2e-16，远远小于0.05，即该因素的影响比较显著，因此就需要进一步利用多重比较方法(如LSD法和Tukey HSD法)来比较该因素不同水平的影响，主要进行判断杨树各个指标实验组和对照组的显著性。

表1杨树方差分析表

由杨树Tukey HSD检验结果(表2)可知，苗高实验组和对照组的均值是21.2156和19.9725，显著性均标记为a；地径实验组和对照组的均值是1.4305和1.3380，显著性均标记为b；生物量实验组和对照组的均值是9.8841和8.0535，显著性均标记为c；可以发现苗高、地径和生物量三个指标显著性差异，因此都可以作为抗涝性指标进行分析。

表2杨树HSD检验信息表

为了进一步检验杨树指标苗高、地径和生物量实验组和对照组差异性，对三个指标分别进行方差齐性检验。由杨树指标方差齐性检验Bartlett's球形检验结果(表)可知，苗高和生物量的P值分别是0.2996和0.6767，均大于0.05，因此无法拒绝原假设，可以进行方差分析；地径的K-squared值是5.7797，P值是0.0162，远小于0.05，因此无法进行参数检验，可以选择非参数检验方法。

表3杨树指标方差齐性检验表

当地径不能通过方差齐性检验时，利用Mann-Withney-Wilcoxon非参数检验比较实验组和对照组抗涝性指标是否存在显著性差异。地径的非参数检验的P值为0.43，远远大于显著性水平alpha值0.05，由此可以认为地径在实验组和对照组中第0天和第60天之间的增量无显著性差异，地径这一指标可以舍弃。因此，在建模分析之前利用苗高变化量和生物量变化量加权构建最终的抗涝性指标Pscore。

利用四分位法删除抗涝性指标Pscore部分离群值，可以实现更好的建模效果，其中离群值筛选公式为

Q

其中Q

表4抗涝性指标描述性统计表

抗涝性指标Pscore的上四分位数是0.4863和下四分数是-0.5196，都比Dscore、Mscore和Sscore其他三个指标要小。剔除异常值(表4)后取值范围更合理，比单一的指标更有研究意义。从图2中可以发现，四个指标的中位数分别是0.09、0.02、-0.12和-0.07，Pscore和Sscore之所以中位数是负数，是因为它们偏小的数值较多一点；散点分布(每个散点用一竖线表示，阴影部分是四分位区间)阴影部分比较短，即Pscore比Dscore、Mscore和Sscore数据更加集中；Mscore的取值范围比Pscore大，这是由于Pscore存在较小的离群点，剔除离群点后Pscore数据集中度明显增加，而且更符合正态分布。

4、杨树品种特征筛选

本实施例所做试验衡量了22个杨树品种的48个相关参数分别在光合作用下和荧光作用下的0天、15天、30天、45天和60天的实验组和对照组的前后具体变化。首先基于以往的研究信息进行专家筛选，剔除了淹水试验中最终死亡的杨树品种(如小叶杨)和与抗涝能力不太相关(基于线性相关性)的特征，得到了20个杨树品种的26个特征，其中还将部分比较相似的特征进行重新构造。如H

专家筛选过程中，还需要考虑特征之间的线性相关性，从专家筛选前后的特征相关性热图对比分析可以发现，筛选之后的剩余特征和抗涝性指标Pscore之间线性相关性更强，如特征qN_Fo、CO2S/R、Fm和H2OS/R的相关系数分别是0.46、0.42、-0.55和-0.47，即特征和杨树抗涝指标Pscore之间线性相关性都是比较高的；过程中同时也剔除了部分线性相关性太高的特征，如Fm、Fv和Fm/Fv三个特征相关性比较高，大量存在会造成数据浪费，将对特征筛选非常不利。专家筛选完成后将更有利于提高后续的模型拟合效果。

接着利用递归特征消除法(RFE)做进一步的筛选。对于RFE函数，其主要参数有：监督学习估计器，主要用来拟合函数，和抗涝性指标比较相关的特征拟合效果越高；系数属性或特征重要性属性是指进行特征选择的相关特征重要性的信息，其取值是给特征打分，分值越高表明其参与度越高，该特征和杨树抗涝性越相关；特征选择数量是指要保留的特征数量，默认保留一半，保留的特征多少也会影响之后的建模分析效果；步数一般默认为1，步数为整数时表示每次要删除的特征数量，当步数小于1时表示每次剔除绝对权重最小的特征。递归特征消除法(RFE)使用一个基模型(基模型可选择一种机器学习模型，如随机森林、决策树等)来进行重复多轮训练：RFE选择好基模型后开始第一轮训练，训练结束后输出杨树样本数据每个特征的系数属性或特征重要性属性，并根据其剔除部分绝对权重较低的杨树抗涝性相关特征，再基于剔除后的特征子集进行下一轮训练。

REF筛选流程可以利用Python中Scikit-learn(sklearn)库具体实现，一般特征筛选数量选择一半即可。本实施例最终选择了在荧光作用下20个杨树品种的13个特征的数据信息参与后续的建模分析。从26个特征中筛选出13个特征如表5所示，分别是AHs/Cs、CO

表5RFE筛选后的特征信息表

5、建模分析与评价

本实施例将杨树抗涝性指标Pscore作为被解释变量,将影响杨树抗涝性的重要特征作为解释变量，采用XGboost、lightGBM、GradientBoost、DecisionTree、RandomForest和Adaboost六种机器学习方法进行定量化分析并分别建立起相应的回归模型。

机器学习方法能够很好地提高模型预测性能，在大量训练之后得到最优模型，但由于中间计算过程属于“黑箱”，无法对变量之间的关系进行直观描述。因此模型的可解释性较弱，需要对杨树抗涝性建立的回归模型需要对模型进行评价，本实施例采用平均绝对误差(MAE)、相对均方根误差(RRMSE)、一致性相关系数(CCC)和纳什萨特克利夫效率系数(NS)、全球绩效指标(GPI)几个指标对各回归模型进行评价分析。

其中，，y

为了能够使用一个标准对本实施例中所有回归模型的预测效果进行排序，将上述四个评价指标归一化得到一个总体指标GPI，GPI表示如下：

其中，GPI

通过上述指标选择出来预测效果最佳和最适合进行杨树抗涝性研究的两个回归模型参与后续的重要特征排名分析。

利用上述四个评价指标分别对各个模型在训练集和测试集的建模效果进行评估分析，其中LightGBM和GradientBoost建立的回归模型表现效果最为优异。利用GPI指数对各个模型在测试集方面进行评价，模型评估结果表明，LightGBM和GradientBoost建立的回归模型的GPI指数最高，而且两个模型在训练集的判决系数都超过了0.96。

6、变量重要性分析

通过GPI指数获得的两个模型分别采用特征重要性(Permutation_importance)和平均SHAP值(Mean Shap Values，通过SHAP方法获得)两种方法对上述筛选出的每一个杨树抗涝性特征进行重要性分析，获得排名前10个的杨树抗涝性特征。

由表6可知，基于LightGBM的模型利用特征重要性(Permutation_Importance)方法进行特征排序(从大到小)的结果是Fm、H

表6两种方法对LightGBM模型特征重要性排序对比

虽然两种方法获取的10个特征的重要性排序有所不同，但10个特征都是完全相同的，这表明采用LightGBM建立的模型的可解释性更强，鲁棒性更高，这意味着对于杨树抗涝性的预测能力也比较强。由Mean Shap Values对LightGBM模型特征排序图(图3)可知，Fm的平均SHAP值最高，其SHAP值散点分布在两侧，都距离0比较远，因而SHAP绝对值的和比较高，而且远远高于其他特征的平均SHAP值；H

由表7可知，基于GradientBoost建立的模型利用特征重要性(Permutation_Importance)进行特征排序(从大到小)的结果是：Fm、RH_S、AHs/Cs、VpdL、Cond、H

表7两种方法对GradientBoost模型特征重要性排序对比

从两种方法对GradientBoost模型特征排序图(图5和图6)可知，Fm、RH_S和AHs/Cs三个特征排名前三，特征重要性远大于其他特征重要性之和。在图6中，Fm、RH_S和AHs/Cs三个特征的盒型长度非常明显即数据分布范围较广，而其余特征的数据取值几乎完全集中在0附近。由Mean Shap Values对GradientBoost模型特征排序图可知，前三个特征的SHAP值有正有负，且远离0的点很多，而其余特征的SHAP值几乎都围绕在0附近；对于每个特征的SHAP值取绝对值后计算平均值便得到平均SHAP值，前三名特征的平均SHAP值依然遥遥领先，比其余特征的平均SHAP值高很多。这将意味着，其他特征的重要性可以忽略，仅仅利用前3个特征进行拟合的预测效果并不会比前10个特征的预测效果差太多，这样一来，根据该模型就可以筛选出比较重要的抗涝性特征。

根据LightGBM、GradientBoost两种机器学习方法建立的回归模型的分析结果能对杨树抗涝性进行很好的预测，这不仅可以节省大量时间成本和运算成本，而且对后续的关于杨树抗涝性的试验设计改进有很大帮助。

以上仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本领域的技术人员在本发明所揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。本说明书中未作详细描述的内容属于本领域专业技术人员公知的现有技术。