基于偏振反射模型和随机森林回归的叶片氮含量反演

张子晗,晏 磊,2,刘思远,付 瑜,姜凯文,杨 彬,刘绥华,张飞舟*

1. 北京大学地球与空间科学学院， 遥感与地理信息系统研究所空间信息集成与3S工程应用北京市重点实验室,北京 100871 2. 桂林航空工业学院， 广西高校无人机遥测重点实验室,广西 桂林 541004 3. 湖南大学电气与信息工程学院,湖南 长沙 410082 4. 贵州师范大学地理与环境科学学院,贵州 贵阳 550001

引 言

在科研人员利用高光谱遥感数据进行植被内部重要物质含量反演之后的30多年里， 高光谱遥感技术已广泛应用于研究与植被相关的多种重要的生物化学物质， 例如： 水分、 纤维素、 蛋白质、 叶绿素、 氮元素、 硫元素和磷元素等。 叶片氮元素含量(leaf nitrogen concentration,LNC)与植物光合呼吸作用和其他的相关生物过程密切相关， 其含量仅占叶片总质量的0.2%～6.4%， 但却能极大程度上影响全球气候变化过程， 并在广域碳氮循环过程中起到主导作用， 对于解构生态环境规律有重大科学意义。 氮元素在叶片内部的主要依托为蛋白质和叶绿素， 这两种物质在可见光近红外(400～2 500 nm)波段有独特的吸收特性， 因而可以利用高光谱数据进行氮元素含量反演。

利用高光谱数据进行氮含量反演的算法可以大致分为三类： 第一类方法建立在叶片反射率与叶片氮含量之间统计关系的基础上， 以逐次线性拟合回归(stepwise multiple linear regression,SMLR)和偏最小二乘回归(partial least squares regression,PLSR)为代表[1]， 在高光谱波段数较多时， 这类算法的回归计算比较耗时， 精度相比于后两类方法更低， 但由于算法原理简单明晰， 仍被广泛使用； 第二类算法建立在高光谱植被指数的基础上， 通过计算高光谱植被指数来反演叶片氮含量[2-3]， 这类算法通过提取红边等深层次光谱特征来构建植被指数， 对于高光谱数据的挖掘仍旧不够深入， 但由于计算简单， 常用于在农学相关的应用领域粗略估计叶片氮含量； 第三类算法建立在机器学习的基础上， 以支持向量机(support vector machine,SVM)、 人工神经网络(artificial neural network,ANN)和随机森林(random forest,RF)为代表[4-5]， 这类算法属于遥感领域与计算机领域学科交叉的范畴， 具有极大的精度优势。 不同的算法选择会产生算法固有误差， 影响氮含量反演精度， 为了尽可能减小算法误差， 同时尽可能利用高光谱数据蕴含的信息， 本研究中测试了多种氮含量反演方法， 最终以反演精度为指标， 选择了第三类算法中的RF算法作为本研究中叶片氮含量反演使用的主体框架。

上述三类氮含量反演算法的输入参数皆为高光谱反射率， 但从定量遥感分析的角度来看， 光谱测量得到的反射率包括两个部分的贡献： 一部分来源于入射光在叶片内部的多次散射过程， 另一部分来源于入射光在叶片表面的镜面反射过程[6]。 在多次散射过程中， 入射光大部分透过叶片表面进入叶片内部， 与叶片内部结构中的叶绿素和蛋白质等含氮生物化学物质发生相互作用， 在叶片内部多次散射后， 携带氮元素信息离开叶片， 这部分来自于叶片内部多次散射的漫反射率所携带的光谱信息是氮含量反演密切关注的。 而在镜面反射过程中， 入射光部分在叶表蜡质层直接发生镜面反射， 这部分反射能量不包含与叶片内部叶绿素和蛋白质等含氮生物化学物质相关的信息， 在氮含量反演中属于高光谱反射率自带的数据误差， 如果使用包含镜面反射部分的高光谱反射率直接进行氮含量反演会产生机理误差， 影响氮含量反演精度。 依据菲涅尔原理， 叶片反射光的镜面反射组分是部分偏振的， 而叶片反射光的漫反射组分由于在叶片内部多次散射导致的消光作用是非偏振的。 因此， 通过偏振反射建模方法获得双向偏振分布函数(bidirectional polarization distribution function,BPDF)， 可模拟计算偏振反射率， 进而扣除测量获取的叶片高光谱反射率中与氮含量反演无关的镜面反射组分， 从源头减少机理误差， 提高氮含量反演精度。

综上所述， 影响叶片氮含量反演精度的误差来源有两个： 算法误差和机理误差。 选择RF算法进行叶片氮含量反演算法可以尽可能减小算法误差， 构建BPDF可以去除大部分机理误差， 最终实现基于BPDF和RF的叶片氮含量精确反演。

1 实验部分

1.1 氮含量反演数据

本研究中用于氮含量反演的数据为主数据集。 叶片氮元素反演的研究地点为德国巴伐利亚国家森林公园[7](49°3′19″N,13°12′9″E)。 根据优势树种的类型， 选择使用国际地圈生物圈计划(international geosphere-biosphere programme,IGBP)土地覆盖分类体系， 优势树种与IGBP类别的对应情况见表1。

表1 巴伐利亚国家森林公园内优势树种对应的IGBP类别Table 1 IGBP classes of dominant vegetation inBavarian Forest National Park

由于本研究区域内的优势树种中阔叶树木均为落叶阔叶品种， 针叶树木都为常绿针叶品种， 所以在研究区域中被简单标注为3种类型： 落叶阔叶林(IGBP04)、 常绿针叶林(IGBP01)和针阔混交林(IGBP05)， 其点位空间分布如图1所示。

图1 叶片氮元素反演点位空间分布Fig.1 Spatial distribution of foliar nitrogen content retrieval plot

共有26个点位作为氮元素反演研究区(8片落叶阔叶林， 8片常绿针叶林， 10片针阔混交林)。 其中， 每个采样区在实地占地为30 m×30 m， 点位范围用分层随机抽样的方法选择8～9棵树木采样， 单颗树木至少采集20片叶片， 妥善保存好后在实验室用凯式定氮法测量氮元素含量。 高光谱数据的获取时间为2013年7月22日， 利用HySpex传感器系统在3 000 m高度获取高光谱影像(HySpex传感器系统包括2台高光谱成像光谱仪， 一台覆盖400～1 000 nm的光谱范围， 另一台覆盖1 000～2 500 nm的光谱范围)， 高光谱数据共有418个窄波段， 410～2 495 nm内， 每5 nm一个通道。 氮含量反演使用的原始反射率是经过几何校正和大气校正处理后的26个样区内的418个波段的双向反射因子数据(bidirectional reflectance factor,BRF)， 除了BRF数据之外无人机系统还记录了获取数据时的观测几何(太阳天顶角、 观测天顶角、 相对方位角)。

1.2 偏振反射建模数据

偏振反射建模数据为辅数据集。 本研究通偏振反射建模获取BPDF， 进而估计偏振反射率。 用于偏振反射建模的BRDF-BPDF数据库[8]由搭载在PARASOL卫星上的POLDER传感器获取。 本研究中为了建立合适的BPDF模型所使用的数据来自于POLDER经过校正的产品， 包含了数据获取时的经纬度和观测几何、 地物的IGBP类别以及地物类型同质性、 样区的归一化植被指数(normalized difference vegetation index,NDVI)、 气溶胶情况、 6个波段(490,565,670,765,865和1 020 nm)的反射率和865 nm波段偏振反射率。

该BRDF-BPDF数据库从POLDER的连续7年的数据中挑选出质量最好的1年， 即2008年， 同时， 在POLDER数据中挑选了质量最优秀(尽可能满足对全球范围的覆盖、 数据有效性高、 样区IGBP类别匀质性高)的点位。 最终数据库的时间跨度和点位选择定为： 包括2008年12个月份的数据， 每个月份采取50个均匀分布在全球的点位， 在每个点位上按照不同的观测几何选择高质量POLDER/PARASOL体系的数据。 本研究使用该数据库中2008年7月份欧洲范围内的数据。

1.3 数据一致性说明

使用两组相对独立的数据集， 主数据集用于氮含量反演， 辅数据集用于建立偏振反射模型， 使用主数据集和辅数据集的流程图如图2所示。 从数据获取的时间来看， 考虑到偏振反射建模具有较强的季节周期性特征， 为了与巴伐利亚国家森林公园的氮含量反演高光谱数据的获取时间契合， 在进行偏振反射建模时仅使用BRDF-BPDF数据库数据中与机载数据同月份的数据， 以保证时间一致性。 从数据获取的地理位置来看， 由于氮含量反演数据的研究位置在德国， 同时考虑到在BRDF-BPDF数据库数据中获取的点位数量不能太少， 最终根据BRDF-BPDF数据库数据的经纬度， 筛选出欧洲范围内的数据进行偏振反射建模， 以保证地理位置一致性。 从两组数据集的空间分辨率来看， POLDER数据的分辨率为6.5 km， 无人机高光谱数据的分辨率为1.65 m(400～1 000 nm)和3.3 m(1 000～2 500 nm)， 但由于BRDF-BPDF数据库中收录的POLDER数据匀质性大于75%， 可认为在数据匹配时受空间尺度的影响较小。

图2 数据使用流程图Fig.2 Flowchart of data utilization

1.3 方法

1.3.1 偏振反射建模

研究中获取的高光谱数据为418个窄波段的BRF。 根据菲涅尔定律， 地表反射是部分偏振的， 其中的偏振部分主要为线偏振， 也是由观测几何(太阳天顶角、 观测天顶角、 相对方位角)决定的， 可用BPDF表征其多角度偏振反射特性。 有文献表明， 在可见光近红外波段， BPDF可视为观测几何稳定时的光谱不变量， 在不同波段的偏振反射率为定值[9]。 BPDF的决定因素是地物类型和观测几何， 为了去除地表反射率的偏振部分， 需要利用BRDF-BPDF数据库中的7月份欧洲范围的数据， 建立3个分别适应IGBP01， IGBP04和IGBP05这3种不同地表IGBP类别的BPDF模型以模拟估计偏振反射率(polarization bidirectional reflectance factor,PBRF)， 进而从BRF中去扣除具有部分偏振特性的镜面反射组分——偏振反射率PBRF， 最终获得漫反射反射率(diffuse bidirectional reflectance factor,DBRF)， 以消除机理误差的影响。

目前精度较高的BPDF模型几乎全部为半经验模型， 纯粹的经验模型虽然能够在局部较好拟合偏振反射率， 但由于缺乏物理机理， 模型参数严重依赖训练数据， 且过拟合明显。 半经验模型通过较少的经验参数对包含物理意义的模型进行校正调整， 其结果是更具有普适性。 本研究选择了5个应用最广、 精度最高的BPDF模型用于模拟IGBP01， IGBP04和IGBP05这3种不同地表IGBP类别下的偏振反射率PBRF， 分别为： Nadal模型[10]、 Waquet模型[11]、 Maignan模型[12]、 Litvinov模型[13]和Diner模型[14]。

为了使公式简洁明了， 在BPDF模型的公式表达中使用一些简记。 其中， 太阳光方向与天顶方向(其夹角称为太阳天顶角θS， 观测方向与天顶方向的夹角称为观测天顶角θV， 太阳光方向在地表的投影与观测方向在地表的投影之间的夹角称为相对方位角φ， 镜面反射过程中的法线与天顶方向的夹角称为半天顶角θh， 入射方向与反射方向之间的夹角的一半为入射角αI， 折射角为αT， 入射方向与反射方向之间方向的变化为方向散射角γ， 折射率为N。 部分角度的余弦值简记为

μS=cosθS

(1)

μV=cosθV

(2)

μI=cosαI

(3)

μT=cosαT

(4)

由菲涅尔公式Fp(γ,N)得到的变量之间的相关关系见式(5)

5种模型的偏振反射建模公式表达见表2。

表2 5种偏振反射模型的公式表达Table 2 Formulas of 5 BPDF models

利用辅数据BRDF-BPDF数据库中的数据分别建立IGBP01，IGBP04和IGBP05这3种地表覆盖类型下的5种BPDF模型， 分别拟合得到自由参数。 通过十折交叉验证的数据分组策略， 以均方根误差(rootmeansquareerror,RMSE)和可决系数(coefficientofdetermination,RSQ)为指标， 挑选3种地表覆盖类型下的最优BPDF模型来模拟偏振反射率PBRF。

1.3.2 氮含量反演

氮含量反演研究基于德国巴伐利亚国家森林公园的实测氮含量主数据、 无人机高光谱主数据和偏振反射建模获得模拟仿真结果， 研究主要包括三个步骤： 首先， 通过BPDF模拟计算偏振反射率PBRF， 进而计算漫反射率DBRF， 获取氮含量反演的算法输入对照组(原始反射率BRF和去偏漫反射率DBRF)； 然后， 设置RF算法参数， 利用整体散射系数和去偏后的漫反射散射系数作为算法输入对照组进行叶片氮含量反演； 最后， 横向对比RF算法和其他算法的表现， 并对偏振反射率在氮含量反演中的影响进行定量评估。

1.3.3 精度评估

偏振反射建模使用POLDER/PARASOL体系的数据库中数据量较大， 采用十折交叉验证的思路进行分组后验证精度， 重复10次十折交叉验证， 取均方根误差RMSE和可决系数RSQ的均值作为BPDF模型精度的度量。 氮含量反演使用26个点位的高光谱数据受到数据量的限制， 采用留一法交叉验证的思路进行数据分组后验证精度， 同样以均方根误差RMSE和可决系数RSQ作为叶片氮含量反演精度的度量， 见式(6)和式(7)。

2 结果与讨论

2.1 BPDF模型选择结果与分析

首先， 按照氮含量反演研究时间、 区域进行辅数据集BRDF-BPDF数据库的筛选， 选择了POLDER数据库中2008年7月采集的欧洲范围内的数据； 然后， 利用十折交叉验证的数据划分思路将数据库中的数据随机划分为训练集和验证集， 利用训练集数据拟合获得5个BPDF模型的自由参数； 获得参数组后利用中位数方法求得单样区最优模型参数组， 根据该模型参数组预测偏振反射率PBRF； 将预测得到的偏振反射率与包含于数据库中的实测865 nm偏振反射率数据进行对比， 最后利用均方根误差RMSE和可决系数RSQ来评价BPDF模型的效果， 5种模型的评价结果见表3。

表3 IGBP01/04/05地表覆盖类型下5种BPDF模型精度评价Table 3 Accuracy assessment of 5 BPDF models on IGBP01/04/05

最终， 以可决系数RSQ为主要评判指标， 均方根误差RMSE为辅助评判指标， 在RSQ尽可能大的情况下选择RMSE尽量小的BPDF模型。 从整体来看， 在同一地表覆盖类型下不同BPDF模型的模拟精度差异不大(保留3位小数后， 均方根误差RMSE差异极小， 可决系数RSQ差异相对明显)， 而在不同地表覆盖类型下， IGBP04地表覆盖类型的BPDF模型表现显著优于IGBP01和IGBP05， 说明阔叶树木的偏振效应相较于针叶树木而言更易模拟， 推测可能与起偏介质的面积大小有关。 模型选择结果为： 在IGBP01地表覆盖类型下最优BPDF模型为Nadal模型， 在IGBP04地表覆盖类型下最优BPDF模型为Litvinov模型， 在IGBP05地表覆盖类型下最优BPDF模型为Maignan模型。 由于不同地表覆盖类型的起偏特性不同， 最优BPDF模型并不统一， 分别利用选出的3种BPDF模型模拟对应的3种地表覆盖类型的偏振反射率。

2.2 RF模型参数设置与叶片氮反演精度横向对比分析

在偏振反射建模的基础上利用原始反射率BRF和漫反射率DBRF作为RF算法输入进行叶片氮含量反演。 RF的核心参数为决策树数量， 本研究分析比对了去偏前后决策树数量为5， 10， 50和100的RF算法的表现， 见表4。 研究结果说明， 在决策树数量一定时， 使用去偏数据进行叶片氮含量反演普遍能获得精度提升且在RF决策树数量设定为10时获得最优氮含量预测模型。 利用RF算法进行氮含量反演能获取的最高精度为： 可决系数达到0.803， 均方根误差达到0.252。

表4 不同决策树数量的RF算法精度评价Table 4 RF algorithm accuracy assessment ofdifferent estimator numbers

2.3 偏振反射率占比统计分析

在26个采样点位， 偏振反射率PBRF在整体反射率BRF中的平均占比为3.580%～11.649%， 各个点位内的偏振反射率占比箱形图见图3， 箱形图形象地展示了偏振反射率PBRF占比的算数均值、 最小值、 下四分位数、 中位数、 上四分位数和最大值。 该分布表明偏振反射率PBRF在整体反射率BRF中平均占比为3.350%， 最高占比16.519%， 说明偏振反射率的影响是不可忽略的， 机理误差足以影响反演精度， 在高光谱叶片氮含量反演中有去除偏振反射率PBRF的必要性。

图3 偏振反射率在整体反射率中占比的分布图中箱体内的×表示算数均值， 箱体内的黑色横线表示中位数Fig.3 Ratio of PBRF to BRFIn this figure,the × in the box represents the mean value,and the transverse line in the box represents the median value

2.4 叶片氮反演精度横向对比分析

横向对比了PLSR算法、 主成分回归(principal component regression,PCR)算法、 支持向量回归(support vector regression,SVR)算法、 K-近邻(K-nearest neighbor,KNN)回归算法和RF回归算法在氮含量反演中的表现， 利用留一法进行交叉验证， 精度评价结果见表5。 其中， 由于偏最小二乘算法应用最为广泛[1]， 将其作为基线算法进行整体精度评估。 研究结果说明， 在利用上述方法进行氮含量反演的过程中， 去除偏振后叶片氮含量反演的精度普遍都有提升(KNN算法在选择的最邻近元数目K足够大时， 最邻近元的选择在去除偏振前后相同， 因而模型预测结果相同， 无精度提升)。 以均方根误差RMSE为核心指标评估提升幅度， 平均提升幅度为4.244%， 其中RF算法去除偏振后叶片氮含量反演精度提升幅度最大， 为13.103%， 且精度最优(达到了上述5种反演算法中的最高可决系数0.803和最低均方根误差0.252)， 充分说明了RF算法利用高光谱数据进行叶片氮含量反演的潜力。 相较于作为基线算法的PLSR(未进行光谱去偏)， 通过偏振反射建模去除大部分光谱机理误差并通过选择随机森林模型减小算法误差后， 叶片氮含量反演精度整体提高32.440%， 证明了去除机理误差和减少算法误差的必要性， 实现了高精度叶片氮含量反演。

表5 5种叶片氮含量反演算法精度评价Table 5 Accuracy assessment of 5 LNC retrieval algorithms

3 结 论

研究旨在提高叶片氮含量反演精度， 从机理误差和算法误差两个误差源入手， 通过偏振反射率建模去除了大部分机理误差， 通过选择精度最优的随机森林算法并调整参数尽可能降低了算法误差， 最终实现了叶片氮含量精确反演， 相比于基线方法整体精度提升为32.440%(可决系数RSQ达到0.803， 均方根误差RMSE达到0.252)。 同时， 在讨论中定量评估了偏振反射率PBRF在整体反射率BRF中的占比分布并对比了去除偏振前后多种氮含量反演算法的精度变化， 体现了随机森林算法在高光谱信息挖掘中的优势， 也证明了去除不携带氮含量信息的偏振反射率的必要性。

致谢：本研究使用的德国巴伐利亚国家森林公园叶片氮含量实测数据和无人机高光谱数据由威斯康星大学麦迪逊分校的王智慧博士提供， 在此表示衷心的感谢！