郭澎涛,朱阿兴,李茂芬,罗 微,杨红竹,茶正早※
(1.中国热带农业科学院橡胶研究所,海口 571101; 2. 农业农村部橡胶树生物学与遗传资源利用重点实验室,海口571101; 3. 海南省热带作物栽培生理学重点实验室-省部共建国家重点实验室培育基地,海口 571101; 4. 中国热带农业科学院土壤肥料研究中心,海口 571101; 5. 南京师范大学地理科学学院,南京 210023; 6. 江苏省地理信息资源开发与利用协同创新中心,南京 210023; 7. 中国科学院地理科学与资源研究所资源与环境信息系统国家重点实验室,北京 100101; 8. Department of Geography, University of Wisconsin-Madison,Madison WI 53706; 9. 中国热带农业科学院科技信息研究所,海口 571101; 10. 海南省热带作物信息技术应用研究重点实验室,海口 571101)
天然橡胶是重要的战略物资和工业原料,其主要来源于橡胶树。磷在橡胶树合成天然橡胶过程中发挥着关键作用,影响天然橡胶的产量和质量。叶片磷含量可以指示橡胶树磷营养状况,因此,获取准确可靠的叶片磷含量是指导管理者合理施用磷肥保障橡胶树健康生长和保持天然橡胶稳产、高产的前提。高光谱技术具有快速、准确估测植物叶片磷含量的潜力。但目前的相关研究主要集中于较小区域,所构建的高光谱估测模型应用范围有限。为增强模型的外推性,一些学者在区域甚至全球尺度上开展植物叶片磷含量高光谱估测研究。但这些研究所构建的高光谱估测模型预测精度并不高,原因是其采集的叶片样本来源于环境异质性较强的区域(气候、成土母质、土壤类型和地形地貌差异明显),受到这些差异很大的环境因素的影响,叶片光谱特征和磷含量发生较大变异,进而导致叶片磷含量与光谱之间关系在区域上产生不稳定性。而这正好与传统全局建模方法所要求的叶片养分含量与光谱之间关系在区域上应当是稳定的假设产生矛盾。为克服这一矛盾,一些学者提出了局部建模方法,该方法假设在变异性较强的大样本集中存在一些局部样本,它们的属性与光谱与之间的关系在空间上是稳定的。若待估测样本与这些局部样本存在相似的目标属性-光谱关系,就可以利用这些局部样本构建高光谱估测模型对待估测样本进行预测。因此,如何选择与待估测样本具有相近或相似目标变量-光谱关系的局部样本就成为利用局部建模方法构建高光谱估测模型的关键。
目前,主要有两类方法用于寻找与待估测样本具有相近或相似目标变量-光谱关系的局部样本。第一类方法为基于样本之间光谱相关性的局部样本搜索方法(Local Sample Searching based on Spectral Correlation,LSS-SC)。该类方法认为被选择的样本与待估测样本之间的光谱相关系数越大,那么它与待估测样本具有相近或相似目标变量-光谱关系的可能性就越大。具体选择过程为先计算待估测样本与光谱数据库中训练样本之间的光谱相关系数,然后设定阈值并选取相关系数高于阈值的训练样本作为局部样本。例如,Shenk等计算每个待估测玉米粒样本与光谱数据库中403个玉米粒训练样本之间的光谱相关系数,并选择大于阈值的训练样本用于构建每一个待估测玉米粒样本的干物质、粗蛋白和酸性纤维含量高光谱估测模型,这些模型的预测精度比传统全局模型提高了6%~13%。该类方法虽然克服了传统全局建模方法存在的缺陷,但容易将与待估测样本具有相似光谱特征但来源于其他环境条件下的样品选中,造成样本误选。
第二类方法是基于样本之间光谱距离的局部样本搜索方法(Local Sample Searching based on Spectral Distance,LSS-SD)。该类方法认为被选择的样本与待估测样本之间的光谱距离越小,那么它与待估测样本具有相近或相似目标变量-光谱关系的可能性就越大。具体选择过程为先计算待估测样本与光谱数据库中训练样本之间的光谱距离,然后设定距离阈值并选择光谱距离低于阈值的训练样本作为局部样本。例如,Ma等计算待估算土壤样点与大光谱数据库中土壤训练样本之间的光谱距离(欧式距离),然后选取光谱距离小于阈值的训练样本作为局部样本构建每一个待估测样本的土壤有机质含量高光谱预测模型,模型预测精度相比传统全局模型有显著提高。该类方法也可用于植物叶片磷含量估测。该方法与LSS-SC相比,虽然在计算光谱相似性时采用了不同指标(光谱距离),但其本质却没有改变,同样仅依赖样本光谱信息,使得该方法依然会从其他环境条件下误选样本。
从上面的分析可以看出,现有的局部样本搜索方法(LSS-SC和LSS-SD)都是以样本之间的光谱相似性为标准来选择局部样本的,只是使用的相似性指标不同。因此,可以将它们称之为基于光谱相似性的局部样本搜索方法。这类方法虽然可以从大样本数据库中搜索到与待估测样本具有相似光谱特征的局部样本,但由于这类方法都是从全部样本中进行局部样本搜索,因此,不可避免地会从其他环境条件下搜索到一些局部样本,而这些样本与待估测样本是不具有相似或相近的目标变量-光谱关系的。利用这些来源于其他环境条件下的局部样本构建待估测样本的高光谱估测模型,毫无疑问会降低模型的预测性能。因此,在利用基于光谱相似性的局部样本搜索方法时如何避免搜索到其他环境条件下的局部样本是需要解决的关键问题。
本研究提出一种先利用环境因子对叶片样本进行类别划分,然后在相同样本类别内部再利用光谱相似性进行局部样本搜索的方法,以解决现有局部样本搜索方法存在的样本误选问题。本文的方法部分将对该方法进行详细陈述,然后在案例研究部分将该方法应用到海南岛北部植胶区,以验证该方法在区域尺度上构建橡胶树叶片磷含量高光谱估测模型的有效性,并将其预测结果分别与利用LSS-SC和LSS-SD搜寻到的局部样本构建的估测模型进行比较。
地理学第三定律指出,地理环境越相似,地理特征(或地理过程)越相近。基于该定律,本研究假设相似环境条件下的橡胶树叶片样本具有相近或相似的磷-光谱关系,且叶片光谱特征越相似,样本之间的磷-光谱关系就越相近。基于这一思想,本研究提出一种基于环境与光谱相似性相结合的局部样本搜索方法(Local Sample Searching based on Environmental Similarity and Spectral Similarity,LSS-ESSS)构建橡胶树叶片磷含量高光谱估测模型。该方法主要包含以下3个关键步骤:1)利用影响橡胶树叶片磷含量的主要环境因子对样本进行类别划分;2)在同一样本类别内利用光谱相似性进行局部样本搜索;3)利用局部样本构建待估测样本的高光谱估测模型。
橡胶树叶片样本类别划分包含2个步骤:1)筛选影响橡胶树叶片磷含量的主要环境因子,2)利用环境因子对叶片样本进行类别划分。
1.2.1 筛选影响叶片磷含量的主要环境因子
影响橡胶树叶片磷含量的环境因子有类别型变量(如成土母质)和连续型变量(如降雨量)。针对不同类型的环境因子采用不同的方式进行变量筛选。类别型变量利用单因素方差分析比较不同类别之间磷含量的差异,若差异达到显著性水平(<0.05),则认为类别型变量对橡胶树叶片磷含量具有显著影响,可被选为影响因子。连续型变量先利用相关性分析选择与叶片磷含量显著相关(<0.05)的变量,然后再利用这些变量拟合与橡胶树叶片磷含量之间的线性回归方程,并依据方差膨胀因子(Variance Inflation Factor,VIF)进行变量共线性诊断。当VIF值大于10时变量之间即存在共线性,需将该变量剔除,并将剩余的变量作为影响橡胶树叶片磷含量的主要环境因子。
1.2.2 利用环境因子对样本进行类别划分
依据所利用的变量类型(类别型或连续型),相应采用不同的样本类别划分方法。若利用的是类别型变量,可直接以不同类别作为划分叶片样本的单元,将位于同一类别内部的叶片样本划分为一类。若利用的是连续性变量,则采用K均值聚类法对橡胶树叶片样本进行类别划分。K均值聚类算法在数据分析软件Matlab 2016a中通过kmeans函数实现。
由于K均值聚类算法不能给出最优的聚类数,因此,本研究采用“手肘法”确定最优聚类数。“手肘法”的含义为所有类别内每个样本与各自质心距离(误差)的平方和(Sum of Square Error,SSE)会随着聚类数的增加而减小,SSE和聚类数之间的关系呈现手肘状,而这个肘部对应的聚类数就是样本的最优聚类数。SSE的计算公式为
式中是聚类数;S代表第个类别;C为S的质心;是属于类别S的样本。
局部样本搜索需要2个步骤完成:1)确定待估测样本所属类别;2)在相同类别内搜索与待估测样本具有相似光谱特征的局部样本。
1.3.1 确定待估测样本所属类别
为实现在相同类别内部进行局部样本搜索,需先将待估测样本分配到相应的样本类别中。若划分类别时利用的是类别型变量,则可依据待估测样本所处的地理位置将其直接划分到相应类别中。若划分类别时利用的是连续型变量,还需要建立待估测样本的类别判别模型。这里以训练集的环境变量聚类结果为基础,结合Matlab 2016a软件中的classify函数建立待估测样本的类别判别模型,并应用该判别模型将待估测样本划分到与其具有相似环境条件的样本类别中。
1.3.2 相同样本类别内部局部样本的搜索
在同一样本类别内利用基于光谱相似性的局部样本搜索方法进行局部样本搜索,即利用LSS-SC或LSS-SD进行局部样本搜索,这两种方法的具体步骤可分别参考Shenk等。
以搜索到的局部样本为建模集,利用偏最小二乘回归(Partial Least Squares Regression,PLSR)以全波谱信息(350~2 500 nm)为输入变量构建待估测样本的高光谱估测模型。PLSR模型中最优潜变量个数通过交叉验证确定,PLSR模型的建立通过R软件中的pls软件包实现。
为验证本研究提出的LSS-ESSS在构建橡胶树叶片磷含量高光谱估测模型方面的有效性,将该方法应用于海南岛北部主要植胶区。在该区域内选择了9个样地(图1)用于橡胶树叶片样品采集。这9个样地的海拔、年平均气温和年平均降雨量的变化范围分别为64~226 m、23.7~24.1 ℃和925~1 773 mm。该区域的土壤类型虽然都为湿润铁铝土,但发育的成土母质却不相同,分别为花岗岩、玄武岩、变质岩和砂页岩,因此,土壤的性质差异很明显。
图1 研究区地理位置及采样点分布Fig.1 Geographical location of study area and distribution of sampling sites
2.2.1 叶片样品
橡胶树叶片磷含量在年内呈现明显的季节变化,4-6月长新叶(抽叶期),叶片磷含量最高;7-9月叶片生长稳定(成熟期),叶片磷含量处于全年中等水平;10-12月叶片衰老(衰老期),叶片磷素逐渐转移到树干和其他组织部位,导致这一时期的叶片磷含量最低。对应于橡胶树叶片生长发育的3个时期,每个时期采集叶片样本1次,每次都在同一地块进行。叶片样品采集时,将每块样地划分为20个区块,每个区块内随机选取5株橡胶树采集叶片,每株橡胶树每次采集树冠下层主侧枝上的没有病斑的顶蓬叶2片。因此,每个区块共采集10片叶子,将这10片叶子作为一个混合样放入一个单独的塑料袋,塑料袋外面记录样品编号、区块地理坐标、种植年限等信息。然后,再把装有混合样的塑料袋放入装有冰块的泡沫箱中。9个样地3次共采集540个叶片混合样。
每次样品采集完毕后迅速送回室内进行光谱测定。光谱测定在暗室内进行,应用的光谱仪为ASD公司的FieldSpec 3可见-近红外光谱仪,波谱范围为350~2 500 nm。在350~1 000 nm范围内,采样间隔和光谱分辨率分别为1.4和3 nm;而在1 000~2 500 nm范围内,采样间隔和光谱分辨率分别为2和10 nm。叶片样品光谱采集时需要先将植被探头通过光纤连接到主机上,植被探头有一内置卤素灯(3.825 V,4.05 W),为光谱测定提供光源。每次测量叶片光谱之前,需要利用叶片夹底部的白板对反射率光谱进行校正,然后将叶片放入叶片夹中进行光谱测定。叶片中部主脉左右两侧区域为测定部位,每个部位测定3次,每个叶片测定6次光谱,一个混合样共采集60次光谱,将这60次光谱进行平均得到混合样品的光谱反射率。
叶片光谱测定完毕后,需要对光谱反射率进行去噪处理。本研究利用Matlab 2016a软件中的butter函数和filtfilt函数进行,这2个函数的计算公式如下:
式中为滤波器的阶数;W为滤波器的截止频率;low代表低通滤波器;和是butter函数返回的滤波系数;spectrum和spectrum分别代表原始光谱和去噪后的光谱。这里构造一个二阶滤波器,即取值2。而W的取值(0<W<1)需要通过试错法确定,W的值越接近0,滤波后的光谱越平滑,反之越接近1,滤波后的光谱越接近原始光谱。本研究尝试了不同的W值(0.9、0.7、0.5、0.3和0.1),发现W为0.9、0.7和0.5时平滑后的光谱依然保留了较多噪音,而W取值0.1时,光谱被过度平滑,一些特征峰消失。W为0.3时叶片光谱的特征峰均能较好体现,同时峰形光滑,说明去噪效果好,因此在这里W取值0.3,滤波后的光谱见图2。构建橡胶树叶片磷含量高光谱估测模型时,以全波谱(350~2 500 nm)作为模型输入变量。
图2 不同时期橡胶树叶片光谱Fig.2 Rubber tree leaf spectrum for different periods
叶片光谱测定完成后,将橡胶树叶片样品放入105 ℃烘箱中杀青30 min,然后降温至70 ℃恒温烘干至恒量,再用研钵磨成粉末过1 mm筛。之后经浓HSO和30%的HO消煮,用钼锑抗比色法测定。化验分析获取的橡胶树叶片磷含量统计结果见图3。
图3 不同时期采集的橡胶树叶片样品磷含量统计结果Fig.3 Statistical results of rubber tree leaf phosphorus concentration for different periods
2.2.2 环境变量数据
有研究表明成土母质和海拔可以显著影响热带森林冠层叶片磷含量。本研究区域中有5种成土母质,分别为玄武岩、花岗岩、变质岩、浅海沉积物和砂页岩,这些成土母质上发育的土壤肥力差异明显。可以预期成土母质可以对橡胶树叶片磷含量产生显著影响。方差分析(图4)证实了这一预想。因此,选择成土母质为影响橡胶树叶片磷含量的主要环境因子。
海拔对森林叶片磷含量的影响主要是通过对气温的影响实现的,即海拔越高气温越低。而在本研究中海拔变化不明显(64~226 m),气温的差异也非常微小(23.7~24.1 ℃)。可以预见在本研究区域中海拔对橡胶树叶片磷含量的影响应该是不显著的,因此海拔未被选择为影响橡胶树叶片磷含量的主要环境因子。
为验证本研究提出的LSS-ESSS的有效性,利用其分别构建每个时期橡胶树叶片磷含量高光谱估测模型。模型构建时,先将每个时期采集的叶片样品随机分割5次,每分割1次就得到1组训练集和验证集,分割5次即得到5组训练集和验证集。样本随机分割时,以每个成土母质为单元以确保每个成土母质中都有训练样本和验证样本。依据每个成土母质中叶片样品数目按比例确定验证集中的样本数。其中,玄武岩有60个叶片样品,从中选取8个样品作为验证样本;花岗岩和变质岩都有40个样品,从各自当中选择4个样品作为验证样本;浅海沉积物和砂页岩都有20个叶片样品,分别从中选取2个样品作为验证样本。这样每个时期的每个验证集有20个样本,每个训练集有160个样本。将LSS-SD和LSS-SC也应用于这些样本中构建橡胶树叶片磷含量高光谱估测模型,并将模型预测精度与本研究提出的LSS-ESSS进行比较。
图4 不同成土母质橡胶树叶片磷含量比较Fig.4 Comparison of rubber tree leaf phosphorus concentration between different parent materials
模型的预测精度通过决定系数(Coefficient of Determination,)、均方根误差(Root Mean Squared Error,RMSE)和测定值标准偏差与RMSE的比值(Ratio of Prediction Deviation,RPD)来衡量。其中,越接近1、RMSE越接近0,RPD的值越大,模型的预测精度越高。上述指标的计算公式如下所示:
为检验本研究提出的LSS-ESSS方法与现有的LSS-SC和LSS-SD构建的橡胶树叶片磷含量高光谱估测模型的预测精度是否有显著差异,利用方差分析对不同模型之间的RMSE和RPD进行比较。
表1列出了利用不同局部样本搜索方法构建的橡胶树叶片磷含量高光谱估测模型在抽叶期的预测精度。可以看出,利用LSS-ESSS构建的模型在第1次、第2次和第3次样本随机分割中的预测精度要高于LSS-SC和LSS-SD,但在第4次和第5次低于LSS-SC和LSS-SD。整体上看,抽叶期利用LSS-ESSS构建的模型预测精度(RMSE=(0.031±0.003)%,RPD=2.027±0.332,=0.697±0.086)要高于LSS-SC(RMSE=(0.034±0.002)%,RPD=1.719±0.158,=0.656±0.053),但差异不显著(>0.05);同时,抽叶期利用LSS-ESSS构建的模型预测精度(RMSE=(0.030±0.004)%,RPD=2.102±0.354,=0.712±0.088)也要高于LSS-SD(RMSE=(0.034±0.002)%,RPD=1.702±0.119,=0.680±0.041),但差异同样不显著(>0.05)。
表1 抽叶期LSS-ESSS模型与其他模型预测精度比较Table 1 Comparison of prediction accuracies between LSS-ESSS and the other models for the period of leaf germination
表2列出了利用不同局部样本搜索方法构建的橡胶树叶片磷含量高光谱估测模型在成熟期的预测精度。可以看出,利用LSS-ESSS构建的模型在5次样本随机分割中的预测精度都要高于LSS-SC和LSS-SD,且在整体上比较(LSS-ESSS的RMSE、RPD和分别为(0.030±0.002)%、2.052±0.196和0.751±0.038,LSS-SC的RMSE、RPD和分别为(0.042±0.002)%、1.491±0.112和0.548±0.052;LSS-ESSS的RMSE、RPD和分别为(0.029±0.003)%、2.202±0.264和0.778±0.043,LSS-SD的RMSE、RPD和分别为(0.042±0.003)%、1.496±0.132和0.559±0.063),差异都达到了<0.05的显著性水平。
表2 成熟期LSS-ESSS模型与其他模型预测精度比较Table 2 Comparison of prediction accuracies between LSS-ESSS and the other models for the period of leaf maturity
表3列出了利用不同局部样本搜索方法构建的橡胶树叶片磷含量高光谱估测模型在衰老期的预测精度。
表3 衰老期LSS-ESSS模型与其他模型预测精度比较Table 3 Comparison of prediction accuracies between LSS-ESSS and the other models for the period of leaf senescence
从表3可以看出,利用LSS-ESSS构建的模型在5次样本随机分割中的预测精度都要高于LSS-SC和LSS-SD,且在整体上比较(LSS-ESSS的RMSE、RPD和分别为(0.026±0.002)%、1.995±0.086和0.760±0.021,LSS-SC的RMSE、RPD和分别为(0.034±0.003)%、1.536±0.120和0.569±0.075;LSS-ESSS的RMSE、RPD和分别为(0.024±0.003)%、2.229±0.143和0.815±0.021,LSS-SD的RMSE、RPD和分别为(0.035±0.003)%、1.492±0.083和0.538±0.053),差异都达到了<0.05的显著性水平。
现有的基于光谱相似性的局部样本搜索方法(LSS-SC和LSS-SD)是从全部样本中搜索与待估测样本具有相似光谱特征的局部样本。这样做的问题是会把其他环境条件下的(如不同母质)样本选作局部样本,这是因为分布在不同环境条件下的样本也可能具有相似的光谱特征。如图5所示,图片左边红圈表示的区域中密集甚至是重叠分布着来自浅海沉积物、砂页岩、变质岩、玄武岩和花岗岩的样本,它们虽然来自不同的成土母质,却具有相似的光谱特征,可称这些样本为混合样本。而来自不同环境条件下的叶片样本往往具有不同的叶片磷-光谱关系,如Asner等在研究亚马逊低地到安第斯山脉这一广阔区域森林叶片光谱变异特征时发现,海拔最高处森林叶片近红外光谱反射率最高,而短波红外反射率最低,此处的氮磷比值最低;但在海拔较低处,森林叶片光谱反射率和氮磷比值变化趋势却与海拔最高处相反。因此,这些混合样本和真正的纯净样本混杂在一起肯定会扭曲叶片磷与光谱之间的关系,导致模型预测性能降低。为避免这种情况的发生,本研究提出基于环境与光谱相似性相结合的局部样本搜索方法,该方法先利用影响叶片磷含量的主要环境因子(本研究中为成土母质)对叶片样本进行类别划分,把来源于不同成土母质的样本区别开来。这时,再从与待估测样本具有相同类别的样本中搜索局部样本,就可以保证搜索到的局部样本与待估测样本具有相同的类别,避免搜索到来自其他成土母质的混合样本。这样就可以保证待估测样本与局部样本具有相近的叶片养分-光谱关系,进而可增强高光谱估测模型的预测精度。
图5 不同时期采集的叶片样本光谱变量主成分分析Fig.5 Principal component analysis for spectral variables of leaf samples collected from different periods
LSS-ESSS适用于大尺度或者环境条件较为复杂的区域。在这些区域,由于环境因子具有较强的异质性,而环境因子又不可避免地会影响到目标变量和光谱特征的变化,最终会导致目标变量和光谱特征在地理空间上也具有较强的变异性,这种变异性会引起目标变量与光谱之间的关系在空间上出现不稳定性。这种关系的不稳定性与传统的全局方法关于目标变量与光谱之间关系在空间上是稳定的假设产生了矛盾,导致传统的全局方法在这种情况下失效。而LSS-ESSS正是针对这一情况发展而来的。
同时,在运用该方法的时候还需要有大量样本的支持。因为该方法需要先将全部样本划分为若干个具有相似环境条件的类别,然后再在与待估测样本具有相同类别的内部搜索局部样本。因此,为了使得搜索出来的局部样本能够充分捕捉到目标变量与光谱之间的关系,就必须要求每一类别内有充足的样本可供选择。所以该方法也特别适合于具有较强数据变异性的大光谱数据库。
针对现有基于光谱相似性的局部样本搜索方法在局部样本搜索时容易将与待估测样本不属于同一环境条件的样本误选的问题,本研究提出基于环境与光谱相似性相结合的局部样本搜索方法,并将该方法应用在案例研究中,得出以下主要结论:
1)来源于不同环境条件下的一些叶片样本呈现出相似的光谱特征,但这些样本不具有相似的叶片磷-光谱关系。
2)本研究提出的LSS-ESSS方法考虑到橡胶树叶片样本所处环境的差异,因此在局部样本搜索时可避免选择到与待估测样本光谱特征相似但环境条件不一致的样本,因而可显著提高局部样本的“纯度”和增强橡胶树叶片磷含量高光谱估测模型预测精度(抽叶期:LSS-ESSS和LSS-SC的 RMSE 分 别 为(0.031±0.003)%和(0.034±0.002)%,LSS-ESSS和LSS-SD的RMSE分别为(0.030±0.004)%和(0.034±0.002)%;成熟期:LSS-ESSS和LSS-SC的RMSE分别为(0.030±0.002)%和(0.042±0.002)%,LSS-ESSS和LSS-SD的RMSE分别为(0.029±0.003)%和(0.042±0.003)%;衰老期:LSS-ESSS和LSS-SC的RMSE分别为(0.026±0.002)%和(0.034±0.003)%,LSS-ESSS和LSS-SD的RMSE分别为(0.024±0.003)%和(0.035±0.003)%)。
3)LSS-ESSS适用于大尺度或环境条件较为复杂的区域,但该方法在应用时需要有大量样本的支撑。