基于主成分分析和随机森林的橡胶小苗叶片磷含量高光谱预测
2016-01-15郭澎涛苏艺茶正早
郭澎涛+苏艺+茶正早
摘 要 为指导橡胶苗施肥管理，通过设置不同磷素水平的橡胶小苗砂培试验，利用高光谱仪测得不同磷处理水平下的橡胶苗叶片光谱反射率，进而运用主成分分析法对原始光谱反射率数据进行压缩及提取主要信息，最后以提取的主要信息为输入变量，结合多重线性回归、偏最小二乘回归和随机森林模型对橡胶苗叶片磷含量进行预测。结果表明：不同磷处理水平下橡胶苗叶片磷含量差异显著；在可见光波谱范围401～669 nm内，光谱反射率随施磷量的增加而增加；不论在训练集还是在验证集，随机森林模型的预测精度最高，训练集和验证集中预测值和实测值之间的相关系数r分别为0.985 0和0.988 4，均方根误差RMSE分别为0.016 5和0.018 2，平均相对误差MRE分别为5.74%和5.99%，模型性能指数RPD分别为3.83和4.01，证明高光谱技术可以快速、准确地诊断橡胶苗叶片磷含量。
关键词 高光谱遥感技术 ；叶片 ；磷素 ；营养诊断
分类号 S794.1
Prediction of Rubber Seedling Leaf Phosphorus Contents Based
on Hyperspectral Technology and Random Forest Approach
GUO Pengtao SU Yi CHA Zhengzao
（Rubber Research Institute，CATAS， Danzhou， Hainan 571737）
Abstract The aim of this study was to test the utility of hyperspectral technology in combination with random forest （RF） approach and principal component analysis （PCA） to estimate rubber seedlings leaf phosphorus contents. A sand culture experiment was carried out to grow rubber seedlings. These rubber seedlings were cultivated with Hoagland's nutrient solutions set at five levels of phosphorus concentration. Leaves of rubber seedlings were sampled at 85， 100， 115， 133， 148 and 163 days after the culture， respectively. A total of 150 samples were collected. At each sampling date， collected leaves were sent to laboratory as soon as possible， and then their leaf hyperspectral reflectance was measured by ASD FieldSpec 3 spectrometer. Phosphorus contents of the corresponding leaves were also analyzed using the conventional chemical analysis method. A second order low-pass digital Butterworth filter with normalized cutoff frequency 0.5 was used to the original spectra to filter out the noise information. Next， PCA was applied to the denoising leaf hyperspectral reflectance to extract the principal components （PCs）. Then，these selected PCs were used as input variables， and multiple linear regression （MLR）， partial least-squares regression （PLSR）， as well as RF model were employed to estimate rubber seedlings leaf phosphorus contents， respectively. Results indicated that rubber seedling leaf phosphorus contents at various treatments were significantly different at p<0.01 level. Within the range 401～669 nm， rubber seedling leaf reflectance increased with increasing leaf phosphorus contents. RF model got the best prediction results. Correlation coefficients （r） between predicted leaf phosphorus contents and measured leaf phosphorus contents were 0.9850 and 0.9884 for train dataset and test dataset， respectively； values of root mean squared errors （RMSE） were 0.0165 and 0.0182 for train dataset and test dataset， respectively； values of mean relative error （MRE） were 5.74% and 5.99% for train dataset and test dataset， respectively， and values of ratio of performance to deviation （RPD） were 3.83 and 4.01 for train dataset and test dataset， respectively. The results verified that hyperspectral remote sensing could be used to rapidly， and accurately predict the leaf phosphorus contents for rubber seedlings.
Keywords Hyperspectral remote sensing technology ； leaf ； phosphorus ； nutrient diagnosis
磷是橡胶树生长过程中必不可少的主要营养元素之一，参与橡胶树体内光合作用、呼吸作用、代谢物质的运输和转运等诸多重要代谢活动。此外，天然橡胶的合成也必须有磷的参与才能完成。缺磷时，橡胶树代谢受阻，茎干增粗变缓，胶乳早凝，严重影响橡胶树的生长和产胶[1]。
在橡胶树营养管理方面，主要通过叶片营养诊断指导橡胶树施肥[2]。传统的橡胶树叶片营养诊断主要依赖实验室化验分析，该方法具有重现性好、准确度高等优点。但其缺点也是显而易见的，即花费大、时效性不强，不能快速、实时指导大面积橡胶树施肥。为解决此问题，亟需发展一种耗费少且快速、准确的橡胶树叶片营养诊断方法。
高光谱遥感技术具有这方面的潜力。国内外相关学者已在利用高光谱遥感技术快速诊断植物氮素营养方面做了大量研究工作，并已形成一些较为实用的模型[3-7]。然而，关于在植物磷素营养方面的研究还较少，目前主要集中在小麦[8-9]、玉米[10-11]、牧草[12-13]、苹果[14-15]和柑橘[16-17]等少数作物上。刘炜等[9]通过分析小麦叶片原始光谱及其一阶导数对不同磷处理水平的响应特征，确定了叶片磷素敏感波长范围，同时提取其对应的吸收面积，并以之为输入变量，结合径向基函数神经网络模型预测了小麦叶片磷含量，预测精度满足了建模要求。Ramoelo等[12]选择淀粉、氮素和蛋白质的光谱吸收波段以及红边位置、简单比率指数和归一化植被指数，并结合一些环境因子为输入变量，利用偏最小二乘回归预测了南非热带稀树大草原牧草氮、磷含量。最近，岳学军等[17]开展了不同生长期柑橘叶片磷含量的高光谱预测研究，利用流形算法对去噪后的光谱数据进行降维和特征提取，进而建立了基于支持向量机的柑橘叶片磷含量预测模型。结果表明，基于一阶导数的Isomap-SVR建模效果最佳，全生长期校正集和验证集模型决定系数分别为0.943 0和0.894 9。
就已有文献可知，目前还未见基于高光谱遥感技术的橡胶树叶片磷含量预测研究报道。上述有关植物磷素营养高光谱诊断的研究都提出了一些所研究作物的磷素敏感波段，但是不同作物之间敏感波段差异较大，这些敏感波段对指导橡胶树叶片磷含量预测意义不大。此外，以单一波段或植被指数为输入变量会造成高光谱信息的损失。另一方面，若以全部原始光谱数据为输入变量进行建模，由于高光谱数据维数远大于目标变量个数，且高光谱不同波段之间存在相关性，导致构建的模型不能准确地反映光谱变量和目标变量之间的关系。同时，植物生化参数与高光谱数据之间的关系往往是非线性的，传统的基于多重线性回归和偏最小二乘回归的预测方法存在一定的局限性。为解决上述问题，本研究引入主成分分析和随机森林法，对砂培培养的橡胶苗叶片磷含量进行预测，以期为高光谱遥感技术在橡胶苗叶片磷素营养诊断方面的应用提供参考依据。
1 材料与方法
1.1 材料
选用经花药组织培养获得的热研7-33-97，供试品种苗龄为5个月。
1.2 方法
1.2.1 试验设计
试验在中国热带农业科学院橡胶研究所5队基地温室大棚进行。试验采用砂培，先将石英砂过80目筛，然后用去离子水冲洗干净并晾干；接着再用去离子水清洗橡胶苗根系，除去培养基质；最后将除去根部培养基质的橡胶苗放入圆柱形塑料培养容器（高30 cm，直径10 cm）中，并用晾干后的石英砂进行覆盖。
橡胶苗用Hoagland营养液进行培养，设置5个磷处理水平： 0（P0）、7.75 mg/kg（P1）、15.5 mg/kg（P2）、23.25 mg/kg（P3）和31 mg/kg（完全Hoagland营养液，P4）。每个磷处理水平设置5个重复（5株苗为1个重复），5个磷处理水平共有125株橡胶苗。
1.2.2 橡胶苗培养
供试橡胶苗于2014年3月2日定植，然后用完全Hoagland营养液连续培养2个月。培养期间，每株橡胶苗每周浇2次营养液，每次浇200 mL。完全营养液培养完毕后，连续14 d用去离子水浇灌橡胶苗，洗去基质内的营养液；接着，利用不同磷处理水平的营养液对橡胶树进行培养，至2014年10月31日结束。期间，每株橡胶苗每周浇灌2次营养液，每次200 mL。
1.2.3 叶片采集与光谱测定
2014年8月11日进行第一次叶片采集，之后每隔15 d采集一次，至2014年10月26日结束，共采集6次。每次采集时，每株橡胶苗取2片稳定期的叶片，每个重复共取10片叶作为1个样本。叶片采集完毕后迅速送回实验室进行光谱测定。
叶片正面光谱反射率用美国ASD公司光谱仪FieldSpec3进行测定。FieldSpec3光谱仪的光谱范围为350～2 500 nm，光谱分辨率在350～1 050 nm为3 nm，在1 051～2 500 nm为10 nm；采样间隔在350～1 050 nm为1.4 nm，在1 051～2 500 nm为2 nm。光谱测定方法详见陈贻钊等[18]。
1.2.4 叶片磷含量化学分析
光谱数据采集完成后，将橡胶苗叶片放入80℃的烘箱中杀青20 min，然后降温至60℃下恒温烘干至恒量，再用研钵磨至粉末，过0.25 mm筛；之后经浓H2SO4消煮，用钼锑抗比色法测定。
1.2.5 原始光谱去噪
原始光谱中存在噪声信息，在应用光谱进行橡胶苗叶片磷含量预测前，需对原始光谱进行滤波处理，以便去除原始光谱中的噪声信息。原始光谱滤波在Matlab8.0软件中运用butter和filtfilt函数进行，butter和filtfilt函数模型分别如下：
[B，A]=butter（N， Wn，‘low）（1）
式中，N为滤波器的介数，Wn为滤波器的截止频率，low代表低通滤波器，B和A是butter函数返回的滤波器系数。在这里N和Wn分别取值2和0.5。
spectrumF=filtfil（B， A， spectrum）（2）
式中，B和A为式（1）中butter函数返回的滤波器系数，spectrum和spectrumF分别代表原始光谱反射率和去噪后的光谱反射率。
1.2.6 主成分分析
主成分分析是一种基于数据统计特征的多维正交线性变化分析。可以将具有相关性的多个变量压缩为几个少数主成分，这些主成分相互独立，且包含原始多个变量的绝大部分信息[19]。主成分分析是在Matlab 8.0软件中通过princomp函数实现的。
1.2.7 模型构建
每个处理随机抽出24个样品，5个处理共抽出120个样品用于建模，余下30个样品用于验证模型精度。本研究利用随机森林（random forest， RF）[20]构建橡胶苗叶片磷含量预测模型，随机森林由分类与回归树（classification and regression tree， CART）发展而来，但与CART又有显著区别[21]。本研究利用R软件中的RandomForest软件包进行RF模型拟合，在对RF模型拟合前，需要对生成树的数量（ntree）、每个节点处用于分割节点的预测变量数（mtry）和叶片最小数量（nodesize）这3个参数进行设置，本研究中，这3个参数分别为1000，6和5。
为了比较RF模型的预测能力，选择多重线型回归（multiple linear regression， MLR）和偏最小二乘回归（partial least-squares regression， PLSR）作为参照模型，MLR和PLSR模型分别在SPSS 13.0软件和The Unscrambler 9.7软件中构建。
1.2.8 模型精度验证
以随机选取的30个样品为验证集，利用验证集中橡胶苗叶片磷含量实测值与模型预测值计算相关系数（correlation coefficient， r）、均方根误差（root mean squared error， RMSE）和平均相对误差（mean relative error， MRE），其中，r越接近于1、RMSE和MRE越小，表明模型的预测精度越高。同时，又选择模型性能指数[22]（ratio of performance to deviation，RPD）对模型的预测能力进行检验，一般认为，RPD≤1.4，表明模型对模拟的对象没有预测能力；1.4
2 结果与分析
2.1 不同磷处理水平下橡胶苗叶片磷含量
表1列出不同采样时期不同磷处理水平下橡胶苗叶片磷含量。可以看出，不同采样时期，随着施磷量的增加，橡胶苗叶片磷含量都呈逐渐增加的趋势。不同采样时期，处理P0与P1之间差异不显著，但P0与P3和P4之间的差异都达到了p< 0.01的显著性水平；除去2014年9月11日和2014年9月29日2次采集的样品，P0与P2之间的差异也达到了p<0.01的显著性水平。
2.2 不同磷处理水平下橡胶苗叶片光谱反射特征
图1为350～2 500 nm内，不同磷处理水平下橡胶苗叶片去噪后的光谱反射率曲线。可以看出，不同磷处理水平下橡胶苗叶片光谱反射率曲线走势基本一致，即在紫外光波谱范围350～400 nm内，随着波长的增加，橡胶小苗叶片光谱反射率缓慢下降；在可见光波段范围401～669 nm内，橡胶小苗叶片光谱反射率先上升后下降，在553 nm处左右形成反射峰；从670～770 nm，叶片光谱反射率急剧上升至一个较高水平；从771～1 266 nm，叶片光谱反射率保持在较高水平上，表现平稳；之后，从1 267 nm开始，叶片光谱反射率逐级下降，且在1 440、1 923 nm处形成2个较深的波谷，在1 660、2 217 nm处形成2个较为明显的波峰。从图1还可以看出，在可见光波谱范围401～669 nm内，不同磷处理水平下光谱反射率随施磷量的增加而增加，而在近红外波谱范围700～2 500 nm内，光谱反射率与施磷量之间的关系并不明显；此外，在整个波谱范围（350～2 500 nm）内，不同磷处理水平下光谱反射率之间的差异不显著。
2.3 橡胶苗叶片高光谱反射率主成分分析
表2为橡胶苗叶片高光谱反射率的主成分分析结果。可以看出，前20个主成分可以解释原始2151条波长（350～2 500 nm）信息的99.993%，因此，选择这20个主成分替代原始2151条波长反射率作为预测叶片磷含量模型的输入变量。
2.4 橡胶苗叶片磷含量预测模型构建
以选取的前20个主成分为输入变量，分别运用RF、MLR和PLSR构建橡胶小苗叶片磷含量预测模型。其中，RF模型为“黑箱”模型，模型不可见，但可以保存和调用。MLR模型如下：
PMLR=0.216-0.026PC1+0.070PC2-0.034PC3+0.124PC4+0.065PC5-0.157PC6+0.283PC7-0.057PC8-0.103PC9+0.092PC10+0.202PC11+0.245PC12+0.989PC13-0.001PC14+0.843PC15-1.092PC16-2.560PC17-0.585PC18+0.181PC19+0.026PC20 R2adj=0.816（p<0.001）（3）
式中，PMLR为橡胶苗叶片磷含量预测值，PC1、PC2、…、PC20为主成分。
PLSR模型如下：
PPLSR=0.21612+2.555PC1-4.331PC2-0.149PC3+5.489PC4-1.885PC5+3.749PC6-7.277PC7+17.954PC8+2.801PC9-0.204PC10+3.352PC11+11.49PC12-22.371PC13-19.763PC14+12.29PC15+58.504PC16-4.216PC17+29.283PC18+49.404PC19
-20.966PC20 R2adj=0.820（p<0.001）（4）
式中，PPLSR为橡胶苗叶片磷含量预测值，其余变量解释同式（3）。
训练集中各模型预测精度可见图2。可以看出，RF模型预测精度最高，其相关系数r和模型性能指数RPD要显著高于PLSR和MLR模型，而均方根误差RMSE和平均相对误差MRE，则明显低于PLSR和MLR模型。
2.5 预测模型精度验证
图3为不同预测模型在验证集中的预测精度。可以看出，RF模型的预测精度最高，其r值和RPD值最大，RMSE和MRE最小。PLSR和MLR模型的r都大于0.9，且其RPD都大于2，MRE值都低于10%，表明这2个模型也都能较准确地预测橡胶小苗叶片磷含量。
3 讨论与结论
3.1 讨论
本研究发现，在可见光波谱范围401～669 nm内，随着施磷量的增加，橡胶苗叶片光谱反射率也逐渐增加。这一发现与Pacumbaba等[23]的研究结果一致。Pacumbaba等[23]发现，缺磷处理的生菜叶片光谱反射率在可见光波谱范围500～665 nm内明显低于施磷处理。但本研究结果与刘炜等[9]、王磊等[10]的研究结果相反，其中，刘炜等[9]在研究缺磷胁迫下冬小麦光谱特征时，发现在可见光波谱范围400～680 nm内，冬小麦叶片光谱反射率随缺磷程度的加深而呈上升趋势；王磊等[10]在对不同磷处理水平下的春玉米光谱响应进行研究时，同样发现在波谱范围350～700 nm内，春玉米大喇叭口期叶片光谱反射率随着磷营养水平的提高而逐渐下降。上述研究结果表明，在可见光波谱范围内，不同施磷处理下不同植物叶片光谱反射率存在明显差异，但光谱反射率随施磷量的变化规律却不相一致，导致这种规律不一致的原因还有待深入探讨。同时，本研究还发现，在近红外波谱范围760～2 500 nm内，不同磷处理水平下橡胶小苗叶片光谱反射率之间的差异不明显，且光谱反射率与施磷量之间无明显的规律。这一结果与Pacumbaba等[23]、刘炜等[9]的结论不一致，其原因有2方面，第一，橡胶树是多年生高大乔木，与一年生的生菜或小麦相比，其对缺磷的响应可能更为缓慢，短期之内不会表现出明显的症状；其次，在近红外波段范围内，植物叶片光谱特征受到水分的显著影响，缺磷状况下的细微光谱特征很容易被叶片水分所遮蔽。
高光谱数据信息量大，且波段之间存在相关性。如果直接将所有高光谱数据用于构建橡胶苗叶片磷含量预测模型，势必会增加模型构建过程中的计算量。同时，由于波段数量远大于样品数量，且波段之间存在共线性，会导致模型估计失真或难以估计准确。本研究通过主成分分析将原始2151维高光谱数据压缩到20个变量（主成分），这20个变量包含原始光谱数据99.993%的信息，且变量之间是两两正交的，即变量之间不存在相关性，这就很好地解决了利用高光谱数据构建模型过程中存在的数据量大和共线性问题。
本研究利用随机森林结合主成分分析法对橡胶小苗叶片磷含量进行了预测，结果要明显优于多重线性回归和偏最小二乘回归，这一结果与Mutanga等[24]的发现相一致，这是因为随机森林可以拟合变量之间的高阶和非线性关系。相关研究表明，植物生化组分与高光谱变量之间并不只是存在简单的线性关系，还存在复杂的非线性关系[25-26]。传统的多重线性回归和偏最小二乘法只能拟合植物生化组分与高光谱变量之间的线性关系，这就限制了其预测精度的提高。这一结果也表明，随机森林在植物生化组分预测方面具有很大的潜力，在未来的研究中可以得到更为广泛的应用。
本研究结果表明，在砂培条件下，应用随机森林结合主成分分析法可以较为准确地预测橡胶小苗叶片磷含量。
参考文献
[1] 何 康，黄宗道. 热带北缘橡胶树栽培：第一版[M]. 广州：广东科技出版社，1987.
[2] 陆行正，何向东. 橡胶树的营养诊断指导施肥[J]. 热带作物学报，1982，3（1）：27-39.
[3] 陈志强，王 磊，白由路，等. 玉米叶片的光谱响应及其氮素含量预测研究[J]. 光谱学与光谱分析，2013，33（4）：1 066-1 070.
[4] 李金梦，叶旭君，王巧男，等. 高光谱成像技术的柑橘植株叶片含氮量预测模型[J]. 光谱学与光谱分析，2014，34（1）：212-216.
[5] Tian Y C， Yao X， Yang J， et al. Assessing newly developed and published vegetation indices for estimating rice leaf nitrogen concentration with ground-and space-based hyperspectral reflectance [J]. Field Crops Research， 2011， 120（2）： 299-310.
[6] Jia F， Liu G， Liu D， et al. Comparison of different methods for estimating nitrogen concentration in flue-cured tobacco leaves based on hyperspectral reflectance [J]. Field Crops Research， 2013， 150： 108-114.
[7] Rotbart N， Schmilovitch Z， Cohen Y， et al. Estimating olive leaf nitrogen concentration using visible and near-infrared spectral reflectance [J]. Biosystems Engineering， 2013， 114（4）： 426-434.
[8] 任红艳，潘剑君，张佳宝. 不同磷肥水平的小麦冠层多光谱特征研究[J]. 土壤，2005，37（4）： 405-409.
[9] 刘 炜，常庆瑞，郭 曼，等. 冬小麦导数光谱特征提取与缺磷胁迫神经网络诊断[J]. 光谱学与光谱分析，2011，31（4）： 1 092-1 096.
[10] 王 磊，白由路，杨俐苹. 春玉米磷素营养的光谱响应及诊断[J]. 植物营养与肥料学报，2007，13（5）： 802-808.
[11] Osborne S L， Schepers J S， Francis D D， et al. Detection of phosphorus and nitrogen deficiencies in corn using spectral radiance measurements [J]. Agronomy Journal， 2002， 94（6）： 1 215-1 221.
[12] Ramoelo A， Skidmore A K， Cho M A， et al. Non-linear partial least square regression increases the estimation accuracy of grass nitrogen and phosphorus using in situ hyperspectral and environmental data [J]. ISPRS Journal of Photogrammetry and Remote Sensing， 2013， 82： 27-40.
[13] Ramoelo A， Skidmore A K， Schlerf M， et al. Savanna grass nitrogen to phosphorous ration estimation using field spectroscopy and the potential for estimation with imaging spectroscopy [J]. International Journal of Applied Earth Observation and Geoinformation， 2013， 23： 334-343.
[14] 朱西存，赵庚星，董 芳，等. 基于高光谱的苹果花磷素含量监测模型[J]. 应用生态学报，2009，20（10）：2 424-2 430.
[15] 潘 蓓，赵庚星，朱西存，等. 基于高光谱的苹果树冠层磷素状况估测模型研究[J]. 红外，2012，33（6）： 27-31.
[16] 黄双萍，洪添胜，岳学军，等. 基于高光谱的柑橘叶片磷含量估算模型实验[J]. 农业机械学报，2013，44（4）：202-208.
[17] 岳学军，全东平，洪添胜，等. 不同生长期柑橘叶片磷含量的高光谱预测模型[J]. 农业工程学报，2015，31（8）：207-213.
[18] 陈贻钊，林清火，华元刚，等. 橡胶树叶片高光谱特征分析[J]. 热带亚热带植物学报，2010，18（3）：293-297.
[19] 朱小鸽. 多重主成分分析及在地质构造信息提取中的应用[J]. 遥感学报，2000，4（4）：299-303.
[20] Breiman L. Random forests[J]. Machine Learning， 2001， 45（1）： 5-32.
[21] Díaz-Uriarte R， de Andrés S A. Gene selection and classification of microarray data using random forest[J]. BMC Bioinformatics， 2006， 7（1）： 3.
[22] Razakamanarivo R H， Grinand C， Razafindrakoto M A， et al. Mapping organic carbon stocks in eucalyptus plantations of the central highlands of Madagascar： A multiple regression approach [J]. Geoderma， 2011， 162（3/4）： 335-346.
[23] Pacumbaba Jr R O， Beyl C A. Changes in hyperspectral reflectance signatures of lettuce leaves in response to macronutrient deficiencies [J]. Advances in Space Research， 2011， 48（1）： 32-42.
[24] Mutanga O， Skidmore A K. Integrating imaging spectroscopy and neural networks to map grass quality in the Kruger National Park， South Africa [J]. Remote Sensing of Environment， 2004， 90（1）： 104-115.
[25] Kokaly R F， Clark R N. Spectroscopic determination of leaf biochemistry using band-depth analysis of absorption feature and stepwise multiple linear regression[J]. Remote Sensing of Environment， 1999， 67（3）： 267-287.
[26] Curran P J， Dungan J L， Peterson D L. Estimating the foliar biochemical concentration of leaves with reflectance spectrometry： Testing the Kokaly and Clark methodologies [J]. Remote Sensing of Environment， 2001， 76（3）： 349-359.
