陈珠琳 王雪峰

(中国林业科学研究院资源信息研究所 北京 100091)

降香黄檀(Dalbergiaodorifera)是海南省特有的极珍贵用材树种，属于国家二级保护植物(郑坚等,2016;苗灵凤,2017),由于具有极高的用材价值，人们对其进行肆意破坏，造成资源急剧减少(杨菁等,2015)。近些年来，海南省大力推广种植珍贵树种，但由于缺乏先进的管理监测技术，推广的效果并不理想。因此，为提高降香黄檀的生长质量与结香速率，必须提出方便快捷的营养诊断方法。传统的营养诊断法需要实验室化学测定，耗费大量人力物力财力(贾良良等,2009)。近些年来，机器自主营养诊断发展迅速。主要方法有：数字图像处理、手持光谱仪和高光谱成像。从数据的丰富度来讲，高光谱成像技术无疑是最佳的选择，但不适用于小尺度监测；手持光谱仪技术可以获取400～2 500 nm之间的植被光谱反射率，但数据之间的冗余度大，易出现休斯现象；数字图像处理技术通常指普通的可见光图像，但由于波段范围窄且带宽较大而受限。相比之下，多光谱相机波段比普通相机丰富，也可以获取各波段的纹理、形状等信息，可用于地面和航空摄影多种尺度，为植物营养诊断提供一种获取数据的新方法。

过量的金属元素富集在土壤中，给植物的光合作用带来不良影响。某些情况下，植被吸收的金属离子会取代叶绿素分子中间的镁原子，减小了聚光并导致光合作用崩溃。重金属污染使植被光谱特性发生明显改变，因此植物的生理敏感现象能够对金属污染做出快速的生理生态反应(屈永华等,2015)。例如，随着铜锌含量的增加，植物在可见光波段的光谱反射率上升，近红外区间的光谱反射率下降，红边出现“蓝移”(王慧等,2017;陈思宁等,2007)。但由于植物种类不同等原因，重金属胁迫导致的光谱反射率特征变化并没有相对一致的结论(Horleretal.,1980;Cleversetal.,2004)。

铁(Fe)是植物生长发育必不可少的微量元素之一，直接参与机体的生理生化活动(袁程等,2012)。尽管土壤中的铁元素总量不低，但在碱性环境中，Fe2+离子含量低，会造成植物缺铁(周晓今等,2012)；而当土壤呈现出弱酸性时，植物又会因“铁毒”出现生理病变。目前，有关叶片金属元素含量预测的研究存在以下不足：1) 大量的研究均将植物叶片或冠层光谱反射率作为自变量，忽视了对纹理特征的分析。2) 输入变量的相关性对模型精度的影响较大，如何减小自变量之间的冗余也是建模的重点。3) 在评价效果相似的模型时，忽略了与植物营养需求的最佳区间相结合，不利于实际应用。鉴于此，本文以幼龄降香黄檀为试验对象，提取多光谱图像特征，主要分析叶片全铁含量与光谱反射率变化，以及植被指数、各波段图像纹理特征与全铁含量之间的关系，提出合理的叶片全铁含量预测方法，为珍贵树种重金属营养诊断提供参考。

1 材料与方法

1.1 试验地概况

研究区设立在海南省文昌市龙楼镇的岛东林场昌洒作业区(19°43′58″—19°44′58″N，110°57′34″—111°01′54″E)，该林场占地约1.9万hm2，属沿海平原地带，热带海洋性气候，年平均雨量为1 808.8 mm，土壤类型为滨海沙壤土，试验地土壤pH 5.0～6.6，有效氮98.3～114.8 mg·kg-1，有效磷3.38～4.56 mg·kg-1，速效钾69.9～78.2 mg·kg-1，有效铁2.33～4.89 mg·kg-1。

1.2 试验设计

试验采用2年生的降香黄檀实生苗，平均苗高为1.07 m，2017年3月选取健壮且生长一致的幼树移植至岛东林场。定植2周后施加等量的水溶性复合肥进行培养，1个月后进行胁迫处理。将试验地划分为4个正方形小区，间隔为5 m，小区面积为25 m2，每个小区内种植25株降香黄檀。单因素试验设计，铁源为乙二胺二邻羟苯基大乙酸铁钠(EDDHA-FeNa)，采用根部施埋法，梯度分别为0、10、15、20 g每株(对应4个小区)，记为CK(对照组)、F1(低铁)、F2(中铁)、F3(高铁)。

1.3 样本采集及多光谱图像获取

试验于2017年4月17日和7月2日进行铁胁迫试验，9月1日进行取样，由于植物体内养分的运输机制导致不同高度的叶片之间养分含量有一定的差异，所以本研究取植株的上中下层，并兼顾内外层，进行后期图像获取和全铁含量化验。

试验使用美国生产的MicaSense RedEdge3多光谱相机获取叶片图像。该相机有5个波段(B、G、R、RE和NIR)，其中心波段分别为475、560、668、717和840 nm；带宽分别为20、20、10、40和10 nm。将叶片放入暗箱中，光源为对应中心波段的单波段LED灯珠组成的灯排，如图1所示。获取叶片图像前进行白板和黑板校正，各通道的光谱反射率计算公式为：

(1)

式中：IC为校正后的反射率，O代表目标物的颜色均值，B代表黑板的颜色均值，W代表白板的颜色均值。

获取图像后，将叶片按编号装入信封并放入烘箱内，105 ℃下烘干0.5 h进行杀青，然后75 ℃烘干至恒质量。将烘干后的叶片磨成粉状，采用原子吸收分光光度法测定全铁含量(简称TIC)。

1.4 自变量提取

叶片图像使用PhotoShop 5.0进行处理，之后提取不同数据源下的图像参数。其中包括：植被指数(以下简称VI)，如表1所示；多光谱图像的各通道纹理特征均值(texture feature mean value,TFMV)与纹理特征方差(texture feature variance,TFV)，包括均值、方差、同质性、对比度、相异性、熵值、二阶矩、相关性。

图1 摄影方法Fig.1 Screen method

表1 植被指数类型及与降香黄檀叶片全氮含量相关性分析①Tab.1 Vegetation indexes and correlation analysis with total nitrogen content in sandalwood leaves

①：*代表在0.05水平上显著相关；**代表在0.01水平上显著相关。下同。* represent significant correlation at 0.05 level;** represent significant correlation at 0.01 level.The same below.

1.5 变量筛选方法

当数学模型的输入自变量多且不相互独立时，会出现过拟合现象，从而导致建立的模型精度低、建模时间长等问题。因此，在建立模型之前，需要将冗余的变量去除，选择最能反映输入输出关系的自变量参与建模。本研究使用4种变量筛选方法：相关性分析法(correlation analysis，CA)、主成分分析法(principal component analysis,PCA)、平均影响值法(mean impact value,MIV)和遗传算法(genetic algorithm,GA)，并比较分析其特点(李哲等,2018;颜胜科等,2017)。

1.6 PSO-BPNN算法

BPNN(back propagation neural network)是目前应用最广的神经网络技术(贾方方等,2016)，具有自适应性、自学习、自组织和非线性映射的功能，适合于模拟各种错综复杂的非线性关系(朱云芳等,2017)，同时也存在收敛速度慢、易陷入局部收敛、网络的学习和记忆能力不稳定等缺点(马廉洁等,2016)。粒子群优化算法(particle swarm optimization,PSO)具有参数少、易实现、精度高、搜索速度快等优点，对多峰问题、非线性均有较好的全局搜索能力(刘耿耿等,2018)，结合BPNN可以解决其学习速度慢以及局部收敛性不足等缺点。将二者相结合，既具有适应性强、泛化能力好的优点，又能减小陷入过拟合的概率，加强预测能力。

1.7 模型验证

本研究按照2∶1的比例随机抽取拟合数据和验证数据。随机抽取65个作为拟合数据，35个作为验证数据。其中模型的拟合优度评价指标为：

(2)

(3)

独立样本检验的统计指标为：

2 结果与分析

2.1 不同铁水平下降香黄檀植株生长量变化

整个试验周期下降香黄檀在树高、冠幅、地茎生长量的变化，结果如图2所示。与对照组相比，施加铁肥不同程度地促进了植株的生长。从CK到F2变化过程中，树高、冠幅、地茎的增长量是持续上升的。但在F3梯度下，树高增长速度与F2相比明显降低，冠幅的增长量则小于F2，但地茎的增长较为明显，增长量是F2的2倍以上。试验结果说明适当添加铁肥有助于提高生物量，但过量后树高、冠幅生长减缓，地茎增大。

图2 不同铁肥梯度下降香黄檀树高、冠幅、地径生长量Fig.2 Growth of height,crown width and stem of Dalbergia odorifera with different iron fertilizer gradientsCK:对照组；F1:低铁梯度；F2:中铁梯度；F3:高铁梯度。CK:Control group;F1:Low level of iron fertilization gradient;F2:Middle level of fertilization gradient;F3:High level of fertilization gradient.

2.1 光谱反射率随TIC的变化趋势

将降香黄檀叶片5个波段的多光谱反射率连接成折线图，如图3所示。NIR波段差异明显，随着施肥量的增加，光谱反射率呈现上升趋势。可见光波段(B、G、R、RE)变化量相对较小，其大小无法识别，图3中将其进行了放大，如灰色矩形框中所示，各波段中不同胁迫对应的反射率排名不同，例如，与其他胁迫组相比，CK在B波段的反射率最高，在G、RE和NIR波段得到的反射率最低。

表2为图2中各点对应的值。B波段光谱反射率在CK～F2梯度持续下降，但在F3梯度出现上升；G波段则相反，即在CK～F2梯度持续上升，在F3梯度下降；R波段反射率则在CK～F1梯度下降，F1～F2梯度上升，F2～F3基本无变化；RE波段则连续上升。这是由于在合理的铁浓度下，叶片叶绿素含量升高，光合作用加强，对B、R、RE波段的吸收加强，G波段的反射率增强；而当全铁出现过量时，植物叶绿体遭到破坏，叶绿素含量降低，则出现与之相反的结果。近红外波段反射率受细胞结构的影响，随着施铁量的增加，植株的生长更加旺盛，细胞分裂与分化能力增强，叶片增大，细胞间空隙增大，所以近红外波段反射率持续增加。

图3 不同全铁含量下降香黄檀叶片光谱反射特征Fig.3 Spectral reflectance of D.odorifera leaves in different TIC

2.2 TIC与图像参数的关系

表1中的皮尔森系数为由不同波段计算得到的植被指数与TIC之间的Pearson系数，与其他波段相比，由RE和NIR波段计算得到的植被指数最能够反映降香黄檀叶片TIC，除MCARIRE和MTVRE外，其他均与TIC在0.01水平上显著相关；而G和NIR得到的植被指数与TIC的现象关系最不明显，仅有4个指数与TIC在0.05水平上相关。R-NIR组合及B-NIR组合得到的植被指数均有7个与TIC在0.05水平上相关，但B-NIR组合在数值上略占优势。

本文对叶片TIC与各波段TFMV与TFV进行显著性检验，结果见表3和表4。TFMV与TIC在0.01水平上相关的有4个，在0.05水平上相关的有11个；而TFV与TIC在0.01水平上相关的有10个，在0.05水平上相关的有12个，这说明变量的离散程度更能反映出叶片TIC值。由表3可知，从波段角度可以得出，RE波段和NIR波段呈现出的相关性最好，二者均在方差和对比度值方面与TIC呈现0.01水平相关，在相异性方面与TIC呈现0.05水平相关。从特征参数角度分析则得出同样的结论，即方差和对比度值相关性最强，其次为相异性。由表4可知：从波段角度得到，NIR波段呈现出的相关性最好，分别在均值、方差、对比度、相异性方面呈现出0.01水平相关。其次是R和RE波段，G波段中呈现相关的参数最多，但均分布在0.05水平。从特征角度分析得到，均值效果最佳，在5个波段中均与TIC具有相关性，其次为相异性、同质性、方差。

综合表1、表3和表4的结论可以得出，RE和NIR波段是指示降香黄檀叶片TIC的重要指示波段，在植被指数以及纹理特征方面均优于其他波段。

表2 不同铁胁迫下降香黄檀叶片各波段光谱反射率①Tab.2 Spectral reflectance of leaves of D.odorifera under different iron stresses

① CK:对照组；F1:低铁梯度；F2:中铁梯度；F3:高铁梯度。CK:Control group;F1:Low level of iron fertilization gradient;F2:Middle level of fertilization gradient;F3:High level of fertilization gradient.

表3 不同波段下全铁含量与纹理特征均值的显著性检验①Tab.3 Significance test between total iron content and mean texture factors in different bands

①NIR:Near IR 波段；RE:红边波段。下同。NIR:Near IR band;RE denotes red edge band.The same below.

表4 不同波段下全铁含量与纹理特征方差的显著性检验Tab.4 Significance test between total iron content and variance texture factors in different bands

2.3 变量相关性分析及筛选

通过公式及图4A可以看出(0～30分别对应表1中的植被指数，顺序为RVIB、RVIG、RVIR、RVIRE、DVIB…，以此类推)，前15个VIs之间的相关性较大，即变量之间具有较强的线性相关性。而后16个VIs之间的Pearson系数多数集中在0.7以下，变量之间的相关性较小。同理对与TIC在0.01、0.05水平上相关的15个TFMV，以及与TIC在0.01水平上相关的10个TFV分别进行相关性分析，结果如图4(b)-(c)所示。图4B中15个TFMV之间的相关性较大，图像呈现出的颜色较深，大部分区域相关性在0.8以上；而图4C中对角线附近的颜色较深，即说明部分相邻参数之间的相关性较大。

图4 参数相关性分析Fig.4 Correlation analysis of parameters

为了减少数据之间的冗余，本文首选选择出与TIC在0.05和0.01水平上相关的变量，再使用4种变量筛选方法进行二次选择。其中相关性分析法中，将VIs和纹理特征的相关性阈值分别设置为0.7和0.8(即比较相关性大于阈值的变量，选择与TIC的Pearson系数最高的变量)；主成分分析保留前三主成分；MIV算法中选择在原值的基础上加减10%构成新的训练样本；GA算法中染色体长度为30，种群大小为20，最大进化代数设置为100。为避免不同自变量个数对最终预测结果的影响，本研究选择每种方法筛选出的前8位变量，最终结果如表5所示。

表5 变量筛选结果①Tab.5 Texture parameter selection results

①VIs代表植被指数；TFMV代表纹理特征均值；TFV代表纹理特征方差。下同。VIs represents vegetation indexes;TFMV represents texture feature mean value;TFV represents texture feature variance value.The same below.

2.4 PSO-BP神经网络模型构建

对每一组试验构建PSO-BP神经网络模型，过程如下：

2) 模型优化。在BP神经网络中内嵌一个种群规模为20，迭代次数为100的PSO算法。为提高算法的全局收敛能力，在PSO算法进化方程中加入惯性权重因子，设置初始惯性权重为0.9，结束惯性权重为0.1，使算法在迭代初期保持较强的全局搜索能力，在迭代后期能进行更精确的局部开发。

2.5 全铁含量预测结果分析

表6 不同试验下PSO-BP神经网络模型预测结果①Tab.6 Prediction results of PSO-BP neural network model in different tests

①CA：相关性分析法；PCA：主成分分析法；MIV：平均影响值法；GA：遗传算法；PSO：粒子群优化算法。下同。CA：Correlation analysis;PCA：Principal component analysis;MIV：Mean impact value;GA：Genetic algorithm;PSO：Particle swarm optimization.The same below.

图5 CA与GA方法实测预测值Fig.5 Measured and predicted value of CA and GA method

2.6 不同特征组合对模型预测结果的影响

由表1～3可知，植被指数对拟合效果的贡献率比较大，所以该参数是必不可少的自变量。本研究以CA和GA为例，设计4组对比试验进行分析，验证纹理特征方面，TFMV和TFV对模型的精度是否有较大程度的提高，结果如表7所示。仅使用VI建模并不能很好的表明降香黄檀叶片TIC，加入纹理特征可以提高预测精度，并且加入TFV的效果要优于TFMV，并且加入2种纹理特征后，各指标提高的幅度较小，这说明纹理特征的离散程度更能表达全铁对叶片的影响。

表7 不同特征组合下的预测结果Tab.7 Prediction results of different feature combination

3 结论

以海南省岛东林场种植的2年生降香黄檀为研究对象，使用多光谱相机获取叶片图像，并建立了反演叶片全铁含量的神经网络模型，得到以下结论：适当添加铁肥有助于提高降香黄檀生物量，但过量后出现树高、冠幅生长减缓，地茎增大的现象。随着叶片含铁量的增大，可见光波段(B、G、R和RE)的反射率变化不明显，NIR波段反射率持续增大。但除NIR波段外，可见光的RE波段也是指示铁胁迫的有效波段，通过上述2个波段计算得到的植被指数和纹理参数与全铁含量的相关性最大。在变量筛选时，使用相关性分析法得到的结果作为自变量得到的预测模型在150～300 mg·kg-1区间更接近真实值，而其他方法在该区间得到的结果偏低，若用于施肥指导，则会增大铁毒产生的概率。另外，仅使用植被指数得到的模型预测结果较差，加入纹理特征后很明显的提高了拟合优度及预测精度，而且纹理特征方差对模型精度的影响更大，即片纹理的离散程度可以较好的作为预测全铁含量的辅助信息。