不同特征筛选方法和估测模型对天然次生林郁闭度估测的影响

王迪赵颖慧

(1.东北林业大学 林学院, 哈尔滨 150040;2.森林生态系统可持续经营教育部重点实验室(东北林业大学), 哈尔滨 150040)

0 引言

森林作为陆地生物圈的主体,在改善生态环境、维护生态平衡等方面起着关键作用。天然次生林占我国森林总面积的70%[1],在森林经营中占有重要地位[2-3]。在森林调查中,郁闭度是不可或缺的因子[4]。郁闭度的测定方法,主要有目测法、树冠投影法、样线法、样点法、鱼眼照片法和冠层分析法等[5],但上述方法仅能获取以样地为单位的抽样数据,难以观测大范围的郁闭度变化情况。遥感技术具有宏观、大面积、多时相和多波段等特点,根据光学遥感数据可对如光谱信息、植被指数[6]和纹理信息[7]等特征因子进行估算,有效地反映森林植被的覆盖情况,但在森林郁闭度高的地区,即使森林复杂度很高,光谱差异还是很小,从而导致光谱信号饱和的现象[8]。

激光雷达(Light Detection and Ranging, LiDAR)所发射的激光脉冲对树冠的响应不仅和高度有关,还与树冠郁闭度有关[9],有效地弥补了光学遥感手段的不足。近年来,有学者使用LiDAR数据提取各种特征郁闭度估测,尤号田等[10]使用小光斑LiDAR点云数据估测樟子松林郁闭度,提取数量比值和能量比值变量建立郁闭度反演模型。穆喜云等[11]使用LiDAR点云密度变量与实测郁闭度进行线性回归。Korhonen等[12]应用LiDAR数据,使用首回波覆盖度指数(First echo Cover Index, FCI)、最后一次回波覆盖度指数(Last echo Cover Index, LCI)和索尔伯格等提出的 Solberg (Solberg’s Cover Index,SCI)指数估测垂直郁闭度和角郁闭度。随着时间的推移,机载激光雷达(Airborne Laser Scanning, ALS)也逐渐用于郁闭度的估算上,Moser等[13]用直角坐标转换为极坐标得到的ALS合成影像计算的郁闭度与从鱼眼照片中计算的郁闭度进行回归,结果表明使用LiDAR合成影像估测角郁闭度精度高、误差小。Ma等[14]利用ALS、机载影像和卫星影像估测郁闭度,进行精度评价,得出ALS数据获得的郁闭度受算法影响最小的结论。

学者们应用各种模型估测郁闭度,其中参数模型主要有多元逐步回归(Multiple Stepwise Regression, MSR)和偏最小二乘回归(Partial Least Squares Regression, PLSR)[15]等;非参数模型有随机森林回归(Random Forest Regression, RFR)、支持向量回归(Support Vector Regression, SVR)[16-17]和Cubist[5]模型等。徐定[15]分别使用MSR、主成分分析(Principal Component Analysis, PCA)、PLSR和像元二分模型估测郁闭度,结合实测郁闭度比较各模型,发现PLSR的精度最高。张瑞英等[18]使用MSR、RFR和Cubist模型估测郁闭度,结果是Cubist模型反演效果最好。

综上所述,尽管采用LiDAR数据及光学遥感数据估测郁闭度的研究已十分常见,但天然次生林郁闭度较高、森林结构复杂,而且结合光学遥感会使植被指数在密林区的敏感性显著下降[19]。因此,本研究以帽儿山实验林场63块样地为研究区域,根据ALS数据提取点云特征,基于3种变量筛选(Pearson相关性分析、随机森林(Random Forest, RF)和Boruta算法)和3种估测模型(偏最小二乘回归(Partial Least Squares Regression, PLSR)、随机森林回归(Random Forest Regression, RFR)和支持向量回归(Support Vector Regression, SVR))估测郁闭度,以鱼眼照片和样点法2种测定方法对估测模型进行精度评价及双因素方差分析,探讨和分析使用ALS数据结合不同特征筛选方法及估测模型对郁闭度较高的天然次生林估测郁闭度的影响,为大范围估测提供了有力支撑。

1 材料与方法

1.1 研究区域概况

研究区为63块(30m×30m)样地,位于黑龙江省尚志市东北林业大学帽儿山实验林场(45°14′～45°29′N,127°29′～127°44′E),平均海拔300 m,如图1(a)所示。植被是典型的东北天然次生林,主要树种有:白桦(Betulaplatyphylla)、蒙古栎(Quercusmongolica)、水曲柳(Fraxinusmandshurica)、胡桃楸(Juglansmandshurica)、黄菠萝(Phellodendronamurense)、山 杨(Populusdavidiana)、色 木 槭(Acermono)、榆树(Ulmuspumila)以及少量人工针叶树种,还有红松(Pinuskoraiensis)、落叶松(LarixOlgensis)和樟子松(Pinussylvestris)等。

1.2 研究数据获取

1.2.1 ALS数据

由运-12飞机搭载LiDAR传感器(Riegl LMSQ680i)和CCD(charge coupled device)相机传感器(DigiCAM-60),于2015年9月在帽儿山林场采集数据,相对航高1200 m。LiDAR传感器波长1 550 nm,最大频率400 kHz,扫描角度在±30°之间,采样间隔1 ns,垂直精度0.15 m,平均点云密度3.6 点/m2。ALS数据预处理包括:点云裁剪、点云滤波、点云分类,使用克里金插值法将地面点插值生成DEM和点云归一化等。

1.2.2 样地获取及2种郁闭度的测定

在帽儿山CCD影像上选取63块30 m×30 m样地,分别于2018年7月和2019年7月进行郁闭度测量,样地分布情况如图1(b)所示,样地信息见表1。研究区为天然次生林样地,为中高郁闭度,植被的生长与更新状况对郁闭度的影响不大,忽略LiDAR数据与样地数据年份不同所造成的影响。本研究采用鱼眼照片法和样点法测定郁闭度。

鱼眼照片法:使用佳能6 D相机和佳能8 mm焦距的圆形鱼眼镜头,在样地中心点和2条对角线的四分之一位置共拍摄5张照片,拍照时尽量选择阴天无云天气,以避免太阳直射光斑对鱼眼照片产生影响,拍摄时注意剔除模糊、受强光照射影响的照片。裁剪掉照片中天顶角大于60°的部分(图1列举的照片即以天顶角60°为准,裁剪后的结果),以排除下木的影响。使用ArcGIS的ISO聚类非监督分类工具进行分类,采样间隔为10像元×10像元,像元共分5类:天空、林隙、树干、树枝和树叶(由于部分照片的天空颜色不统一,因此将有云的分为7类,没有天空颜色或云干扰的分为6类或5类),输出分类栅格,对比分类结果图和原图,分别统计林冠像元及总像元个数,二者之比为郁闭度,每个样地计算的5张鱼眼照片郁闭度均值作为鱼眼照片郁闭度。

表1 样地信息统计表Tab.1 Statistical table of plot information

图1 帽儿山实验林场位置、63块样地分布图及部分鱼眼照片示例图Fig.1 The location of Maoershan Experimental Forest Farm, distribution of 63 sample plots and parts of fisheye photos

样点法:沿样地对角线一步一抬头,判断样点是否被树冠遮盖,郁闭度为被树冠遮盖的样点数和样点总数的比值。

1.3 特征因子提取

以ALS数据计算首回波的高度和强度为指标,即均 值(Mean)、方 差(Variance, var)、标 准 差(Standard deviation, std)、最小值(Min)、最大值(Max)、中位数(Median)、变异系数(Coefficient of variation, cv)、偏度(Skewness, skew)、峰度(Kurtosis, kurt)、下四分位数(Q1)和上四分位数(Q3)等,共22个特征。另外,计算首回波覆盖指数(公式中用FCI表示)和激光穿透指数[20](Laser Penetration Index, LPI,公式中用LPI表示),公式为

式中:Scanopy为只有一次回波的冠层点云数;Fcanopy为首回波的冠层点云数;SAll为只有一次回波的全部点云数;FAll为首回波的全部点云数;Iground为地面点云总强度;Ivegetation为冠层点云总强度。

1.4 特征筛选方法

特征选择算法可以有效避免维度问题[21]或Hughes效应[22],从而改进估测性能[23],因此本研究采用3种特征变量筛选方法。

1.4.1 相关性分析

相关性分析(Correlation Analysis)是描述2个变量之间关系紧密程度的指标[24]。本研究用SAS 9.3进行相关性分析,计算自变量与郁闭度的Pearson相关系数。

1.4.2 随机森林

随机森林(Random Forest, RF)在特征选择减少冗余数据方面效果很好,本研究在 python 3.7中进行RF筛选,设置2个参数:决策树数量(n_estimators)为1 000,随机状态(random_state)为6,使用平均下降精度(Mean Decrease Accuracy, MDA)来评价特征变量的重要性。

1.4.3 Boruta算法

Boruta算法是基于RF算法思想构建的特征选择算法,主要通过循环方法评价特征变量的重要性,对每个特征值随机混合构造具有随机性的阴影特征[25];以在RF算法的每次迭代中,得到原始特征和阴影特征重要性的Z分数最大值(Max_Shadow)为筛选指标,特征变量Z分数大于Max_shadow时,保留该特征,从而筛选出最优特征集合[26]。

1.5 郁闭度估测

1.5.1 偏最小二乘回归

偏最小二乘回归集多元线性回归分析、PCA和典型相关分析的基本功能为一体,可以解决自变量之间存在的多重共线性问题[27]。

1.5.2 随机森林回归

随机森林回归[5]建立在决策树基础上,通过多次bootstrap抽样获得多个随机样本,并通过这些样本分别建立对应的决策树,从而构成RF。在回归过程中,放回随机抽取样本,因此建立回归树时部分样本不会被选中作为检验样本出现,起到了样本内部交叉验证的作用,2个随机性引入可减少过拟合情况的发生。

1.5.3 支持向量机回归

支持向量机回归[28]的基本思想是将实际问题按照非线性对应关系映射到高维空间,再进行线性回归的处理,最终得到原始空间的非线性回归结果。本研究使用Python3.7的sklearn包来实现,核函数为线性核函数,多项式系数为3,误差项的惩罚参数C=1,损失函数的值为0.1。

1.6 精度评价

本研究仅有63块样地,样地数量少,因此采用留一交叉验证法(Leave-One-Out Cross Validation, LOOCV)[28],使用决定系数(R2)、平均误差(Mean Error, ME)、平均绝对误差(Mean Absolute Error, MAE)和均方根误差(Root Mean Square Error, RMSE)4个指标[5,29]对模型进行精度评价。

1.7 方差分析

本研究目的是探讨特征变量筛选方法和估测模型对天然次生林郁闭度估测的影响,因此在3种特征变量筛选方法(Pearson相关性分析、RF和Boruta算法)和3种估测模型(PLSR、RFR和SVR)作用下,利用SAS9.3进行方差分析(Analysis of variance, ANOVA),多重比较采用LSD法(Least Significance Difference, LSD)[30]。

2 结果与分析

2.1 特征变量筛选结果

2种郁闭度测定方法(鱼眼照片法和样点法)的相关性分析结果见表2,2种方法均保留了15个变量且所保留的变量有很大的一致性,点云的高度和强度与2种测定方法相关性都很高。由表2可以看出,高度变异系数(hcv)与2种测定方法的相关性最高;FCI和LPI与2种方法的相关性均在0.6左右,这与其计算公式(公式(1)和公式(2))有关,FCI为正相关,LPI为负相关;强度变异系数(icv)和强度峰度(ikurt)与鱼眼照片法相关性不高(低于0.3),强度偏度(iskew)和ikurt与样点法的相关性均低于0.3。

表2 鱼眼照片和样点法的Pearson相关性分析结果Tab.2 The results of Pearson correlation analysis based on fisheye photo and sample method

图2为2种郁闭度测定方法(鱼眼照片法和样点法)的RF特征变量重要性排序结果。2种方法均保留了14个特征变量,保留的变量大部分相同,但重要性不同,鱼眼照片的FCI、hcv和hskew重要性最高,保留了hmean、ivar、高度最小值(hmin)和hmedian等特征变量;样点法的hskew、iQ3和hcv重要性最高,保留了高度方差(hvar)、高度标准差(hstd)、imedian和强度最小值(imin)等特征变量。

图2 鱼眼照片和样点法的RF变量筛选结果Fig.2 The feature selection results of RF based on fisheye photo and sample method

图3为2种郁闭度测定方法(鱼眼照片法和样点法)的Boruta算法变量筛选结果。在Boruta算法变量筛选中,Z分数值大于Max_Shadow(Z分数最大值)的变量作为接受变量被保留下来。2种方法经过Boruta算法筛选后保留的特征变量有所区别。鱼眼照片法的LPI、imean、FCI和hQ1作为接受变量被保留下来;样点法的hmean、hcv、FCI和hQ1作为接受变量被保留下来。

图3 鱼眼照片和样点法的Boruta算法变量筛选结果Fig.3 The feature selection results of Boruta algorithm based on fisheye photo and sample method

2.2 郁闭度估测模型精度评价

表3为2种郁闭度测定方法(鱼眼照片法和样点法)的模型精度评价结果。鱼眼照片法中,PLSR精度最高,SVR最低。Boruta算法筛选后的PLSR估测精度最高(R2=0.451 1,RMSE =0.067 5),Pearson相关系数筛选后,SVR估测精度最低(R2=0.303 4,RMSE =0.079 1)。另外,ME和MAE值越小,表明预测值与实测值差距越小,估测效果越好,表中PLSR和RFR的ME值大多数小于0.1,SVR模型存在严重的低估现象。对于样点法,RF筛选后的RFR估测精度最高(R2=0.372 9,RMSE =0.079 2),Pearson相关系数筛选后的SVR估测模型精度最低(R2=0.062 3,RMSE =0.096 9)。2种估测模型(PLSR和RFR)的ME值多数小于0.1。经Boruta算法变量筛选后的SVR模型的ME值为0.519 6,存在着高估现象,另外2种筛选方法的ME分别为-0.315 4和-0.417 7,模型存在着低估现象。样点法的SVR模型拟合情况远不如鱼眼照片法,虽然经Boruta算法变量筛选后SVR模型(R2=0.192 2)要高于其他2种筛选方法(0.062 3和0.106 5),但与2种郁闭度测定方法的其他模型相比,估测效果是最差的。

表3 鱼眼照片法和样点法对估测郁闭度模型精度评价结果Tab.3 The results of accuracy assessment of estimation models based on fisheye photo and sample method

绘制2种方法中估测精度最高的变量筛选方法郁闭度实测值与估测值拟合曲线,图4(a)—图4(c)为鱼眼照片法拟合曲线,图4(d)—图4(f)为样点法拟合曲线。从图4(a)—图4(c)可以看出,鱼眼照片的中低郁闭度区域内,存在明显的低值高估现象;高郁闭度区域内,高估低估情况都有发生,但高值低估现象更为严重。由图4(d)—图4(f)可以看出,样点法经RF变量筛选后RFR和PLSR,在中低郁闭度区域内,存在低值高估现象,但使用SVR,则出现了低估的情况;高郁闭度区域内,高估低估情况都有发生,但高值低估现象更严重,在SVR模型中,数据的散点比较分散,估测效果不佳。

图4 鱼眼照片和样点法中估测精度最高的变量筛选方法郁闭度实测值与估测值拟合曲线图Fig.4 The fitting curve between the measured and estimated values of canopy closure of feature selections method with the highest estimation accuracy in fisheye photo and sample method

2.3 方差分析结果

使用SAS9.3对数据进行方差分析,分析特征筛选方法和估测模型2种因素对郁闭度估测的影响及其交互作用。由表4可以看出,鱼眼照片法的P均大于0.05,变量筛选方法和模型拟合方法对郁闭度估测无显著影响,也不存在交互作用。样点法的变量筛选方法对郁闭度估测无显著影响,模型拟合方法有显著影响,变量筛选方法和模型拟合方法之间不存在交互作用。

由模型均值比较图5看,样点法的SVR模型相比于PLSR和RFR有显著差异。

表4 鱼眼照片法和样点法的方差分析结果Tab.4 The results of ANOVA based on fisheye photo and sample method

图5 样点法的3个模型的均值比较图Fig.5 Comparison of the mean values of the three models of sample method

2.4 郁闭度分布图

以估测效果最优的Boruta算法筛选后的PLSR模型为例,制作帽儿山实验林场郁闭度空间分布图,如图6所示。由图6可以看出,帽儿山实验林场的郁闭度较高,大多数在0.8左右。

图6 帽儿山实验林场郁闭度分布图Fig.6 Canopy closure distribution of Maoershan Experimental Forest Farm

3 讨论

3.1 特征变量筛选

由3种筛选方法保留的特征变量可以看出,FCI、LPI和hcv既是与郁闭度相关性最高的变量,也是重要性较高的变量,FCI(冠层点云数与所有点云数的比值)和LPI(地面点云与所有点云的强度比值)的计算公式(公式(1)和公式(2))与郁闭度的概念最接近,而hcv侧重表达冠层垂直结构的波动程度,这3个特征变量对郁闭度的解释性更强,对模型的贡献最大,而高度均值代表了冠层中上层信息,hQ1等变量同样对模型起到了很好的解释作用。

对Boruta算法来说,尽管保留的特征变量数远少于其他2种方法(Pearson相关性分析和RF),但保留的变量与郁闭度有很强的相关性和重要性。阴影特征的引入,给Boruta变量筛选增添了更大的随机性,既能更好地选出对模型重要性更高、贡献度更大的变量,也起到了降维作用,提高了模型运行的效率,同时很好地保证了模型的精度,这与卢永亮等[25]的研究结果相同,同时本研究还对比了3种变量筛选方法,发现Boruta算法筛选的变量,参与模型拟合时精度最高(R2=0.451 1,RMSE =0.067 5),体现了其在变量筛选上的优势。

3.2 郁闭度估测模型

3种估测模型的精度,大部分情况下PLSR的模型估测精度高于RFR,SVR最低。本研究提取的点云高度和强度统计量之间具有很高的相关性,PLSR模型克服了传统的多元回归模型的缺点,避免了多重共线性的影响,在本研究中拟合效果最好。RFR的估测精度略低于PLSR,可能是用于模型拟合的自变量过多,导致了过拟合情况,也可能是样本数目较少对模型训练结果产生影响。帽儿山林场天然次生林的植被覆盖率很高,低郁闭度的样地数量不多,大部分的样地均为中高郁闭度,数量较多的高郁闭度样地数据对回归模型的拟合产生影响,导致低郁闭度样地的模型拟合效果不佳,误差较大,而SVR比RFR估测精度低,可能是因为核函数在不同的研究领域适用性不同,所选参数对数据拟合也有一定影响,在使用样点法测定郁闭度时,尽管已经使用相关性高以及重要性高的变量,SVR的精度却远不如其他2种模型(虽然经Boruta算法变量筛选后精度最高,R2也仅为0.192 2),说明SVR并不适用于本研究中样点法测量的郁闭度。郁闭度估测时,3种模型都出现低值高估和高值低估的现象,这与张瑞英[18]的研究一致,主要因为天然次生林中低郁闭度的样地数目少,代表性不强,模型主要反映高郁闭度的样地情况。

3.3 不同特征变量和估测方法对估测郁闭度的影响

在2种郁闭度测定方法(鱼眼照片法和样点法)中,大多数情况下,由于Boruta算法筛选后的变量相关性和重要性都很高,所以Boruta算法筛选变量的效果优于RF和Pearson相关系数法,因此,将Boruta算法筛选后的变量引入回归模型时,模型的精度更高。同时,PLSR模型能很好地克服多重共线性,拟合效果更好,所以使用鱼眼照片法计算郁闭度时,经Boruta算法筛选变量后的PLSR模型精度最高(R2=0.451 1,RMSE =0.067 5)。

从2种郁闭度测定方法来看,鱼眼照片法的估测精度大多数情况下都优于样点法。当使用PLSR时,经Pearson相关性分析变量筛选后,鱼眼照片的R2比样点法高16.44%,RMSE降低了14.18%;RF变量筛选后的R2提高17.28%,RMSE降低了14.47%;Boruta算法变量筛选后R2提高了17.8%,RMSE降低了14.88%。当使用RFR时,除了Pearson相关性分析筛选,鱼眼照片法的R2比样点法低4.49%,RMSE降低了7.87%;RF变量筛选后的R2提高8.6%,RMSE降低了11.49%;Boruta变量筛选后的R2提高了7.29%,RMSE降低了10.93%。当使用SVR时,鱼眼照片法的估测精度尽管略低于其他模型,但仍处于较为正常的水平,但使用样点法测定郁闭度时,SVR的估测精度与PLSR和RFR有明显差异,Pearson相关性分析的R2为0.062 3,RF变量筛选后的R2为0.106 5,Boruta算法变量筛选后精度最高,R2也仅达到0.192 2,且2种郁闭度测定方法中,用于模型回归的自变量差异不大,说明本研究中基于鱼眼照片法测定的郁闭度的模型表现更好,在郁闭度估测上比样点法更为精准,这与智献坡[31]的研究结果相同,但与何萍等[32]的研究结果有所不同,可能是样点法测定郁闭度时,仅在样地对角线上记录样点,导致样点数目少,样点法估测效果不好,应适当增加样地内的样点分布。

4 结论

本研究使用ALS数据和2种样地郁闭度测定方法(鱼眼照片法和样点法)探究3种特征筛选方法(Pearson相关性分析、RF和Boruta算法)和3种模型(PLSR、RFR和SVR)对天然次生林郁闭度估测的影响。结果表明:RF和Boruta算法变量筛选后的估测模型精度要高于Pearson相关系数。PLSR克服了传统的多元回归模型的缺点,避免了多重共线性的影响,拟合效果最好,SVR最低。鱼眼照片法比样点法精度高,其中Bortua算法变量筛选后以PLSR估测郁闭度的精度最高(R2=0.451 1,RMSE =0.067 5)。样点法和鱼眼照片法中的变量筛选及鱼眼照片法中的模型对估测郁闭度无显著影响,样点法的模型拟合方法对估测郁闭度有显著影响,SVR模型与PLSR、RFR有显著差异。