基于GF-1数据多尺度遥感特征的森林蓄积量估测研究

黄冰倩,岳彩荣,朱泊东

(1.西南林业大学,昆明 650224；2.贵州省林业调查规划院,贵阳 550003)

随着遥感(RS)、地理信息系统(GIS)、全球导航卫星系统(GNSS)等信息技术在林业资源监测应用日趋成熟,国内外已有一些学者利用“3S”技术结合数学模型对蓄积量进行遥感估测,并尝试了一些新的方法[1-2]。遥感影像早期较多利用Landsat TM/ETM,SPOT5,IKONOS提取遥感特征因子,建立森林蓄积估测模型。遥感特征因子也从波段光谱特征、植被指数特征逐渐增加了地形特征、纹理特征因子[1-2]。Shamsoddini[3]应用WorldView-2影像,对澳大利亚新南威尔士州的松树林分蓄积进行了估测。近年来,我国自主研制的卫星遥感数据(如高分1号,高分2号,资源3号等)逐渐推广应用,且影像时间分辨率、空间分辨率、覆盖区域都能较好地满足研究要求。结果显示，波段的纹理特征比光谱值和衍生光谱值具有更好的反演效果[4-6]。通过随机森林方法对不同的特征变量重要性进行研究,发现纹理特征比光谱特征因子对森林蓄积建模影响大[4-6]。随着引入特征变量的增加,一方面提高了建模精度,另一方面也会引起共线性问题。通过多元逐步回归、随机森林建模的方法,筛选关键变量,提高建模精度,减少数据冗余。

遥感森林蓄积估测模型研究逐渐由参数化方法向非参数方法模型过渡,参数化方法通用性、可解释性较好,但是建模精度比非参数方法低,且建模因子较多时会存在数据冗余等问题[7-11]。与参数化方法相比,非线性森林蓄积估测模型拟合精度较高,无需先验假设,模型构建更便捷。机器学习算法是典型的非参数化模型构建方法,虽然此类算法不能输出具体模型,但也可以进行回归预测[7-11]。目前,通过不同因子组合、对比不同建模方法来研究模型拟合精度较多,从特征提取单元与地面监测尺度匹配性考虑较少。本文基于贵州省观山湖区GF-1遥感数据提取的光谱信息、植被指数及纹理特征,结合研究区实测马尾松样地数据,通过多元逐步回归、随机森林算法构建不同窗口遥感特征的森林蓄积量估测模型,寻求最优拟合效果对应的适宜窗口,进一步深化提高森林蓄积估测的建模精度方法,旨在为高分辨率光学遥感影像数据估测森林蓄积探索一些方法和思路。

1 研究区概况

研究区位于贵州省贵阳市中心城区西北部的观山湖区,地处苗岭山系西部中段,东邻云岩区黔灵镇、南接花溪区久安乡、西靠清镇市麦格苗族布依族乡、北接观山湖区艳山红镇。地理位置 26°33′～26°40′N,106°33′～106°41′E,国土总面积306.95km2。地貌以中山丘陵为主,平均海拔1 280m左右,地势相对平缓。地带性植被属于亚热带常绿阔叶林,植被大多是次生植被和人工植被,主要树种有马尾松、杉木、柳杉、泡桐、青冈、麻栎等。

2 研究方法

2.1 数据来源及处理

影像数据来源于2015年10月观山湖区高分1号(GF-1)遥感影像(图1),全色波段分辨率为2m,多光谱波段分辨率为8m,投影为2000国家大地坐标系(CGCS2000)。该影像为开展二类调查统一采购,已经过辐射定标、大气校正、几何校正、影像拼接、匀光匀色等处理,并通过假彩色合成得到观山湖区范围内高分遥感影像。GF-1光谱对应情况如表1所示。研究区影像完成几何配准后,对全色和多光谱影像进行融合和裁切,输出2m分辨率的GF-1融合影像。在融合影像基础上,生成分辨率为2m的各波段影像(蓝、绿、红、近红外),并以此为基础数据进行后续的光谱特征、植被指数计算、纹理特征提取。

图1 研究区高分1号遥感影像图

表1 高分1号卫星波段光谱信息

2.2 样地布设

外业数据主要采用2016年观山湖区二类调查中马尾松样地实测数据,样地布设在GIS平台上,以1∶10 000遥感影像图为基础,在调查区域内随机确定起点后,以计算出的样地间距(L)为参照,从上至下、从左至右按公里网布设样地,样地规格为半径14.57m的圆形样地,面积667m2。对样地内胸径大于等于5.0cm的样木进行每木检尺,采集其胸径、树高、郁闭度、优势树种等数据,对各树种(组)检尺样木按径级组归类,用各径级组平均树高、各径级组检尺木胸径通过二元立木材积式计算各株样木材积,汇总得到样地蓄积,记载到0.1m3。研究区选取马尾松样地点总样本共180个(其中训练样本120个,验证样本60个)各样本最大值、最小值、平均值情况及样地空间分布情况如表2所示。

表2 研究区样地蓄积量分布统计表

2.3 特征变量设置

2.3.1光谱特征

1)光谱信息。以高分1号影像的band4(近红外)、band3(红光)、band2(绿光)、band1(蓝光)波段作为单波段遥感因子。在ArcGIS中使用“区域分析”中以“以表格显示分区统计”工具提取特征单元对应遥感数据的反射率。

2)植被指数。在ArcGIS中利用字段计算器,进行植被指数计算。

归一化植被指数:NDVI=(B4-B3)/(B4+B3)

比值植被指数:RVI=B4/B3

差值植被指数:DVI=B4-B3

式中:B4为高分1号影像的第4波段；B3为高分1号影像的第3波段。

2.3.2纹理特征

纹理特征主要体现遥感影像的纹理信息,主要包括均值(ME)、协同性(HO)、方差(VVA)、相关性(CO)、二阶矩(SM)、相异性(DI)、熵(EN)、对比度(CT)。本文基于ENVI软件选用二阶概率统计的滤波方法,在3×3,5×5,7×7,9×9,11×11,13×13,15×15,17×17,19×19,21×21等9种窗口下分别提取GF-1各波段数据的纹理特征。

2.4 森林蓄积估测

选用多元逐步回归方法、随机森林算法构建森林蓄积量估测模型。

2.4.1特征变量分析

在多元逐步回归建模过程中,参与拟合的自变量越多,建模精度将会提高,同时对建模无显著影响的因子可能会引起共线性和数据冗余问题,从而影响估测结果。因此,对森林蓄积量和遥感因子进行相关性分析,根据决定系数R2、调整R2、精度检验指标等选取主要变量是线性回归建模的关键。随机森林的特征选择方法属于包裹法,可在训练过程中输出变量的重要性,并对特征变量重要性进行排序,即哪个特征变量对预测类更有用。经多次试验确定决策树数目和节点分裂时变量个数,代入随机森林回归模型,并利用回归模型对验证样本进行预估[12-17]。

2.4.2窗口与样地的匹配

滤波方法提取影像的纹理特征主要受窗口大小、X和Y变换值、选择灰度量化级别等参数的影响。在提取影像纹理特征的过程中,设置X和Y变换值默认值为1和1,灰度量化级别默认值为64,对比设置3×3,5×5,7×7,9×9,11×11,13×13,15×15,17×17,21×21不同窗口。本研究样地规格半径为14.57m,面积667m2的圆形样地,使用GF-1单波段遥感影像单个像元面积为2m×2m,以13×13的移动窗口取平均值,代表676 m2遥感特征值,与实测的圆形样地面积尺度上可较好匹配。

2.4.3模型的构建与拟合效果

通过特征变量筛选、最佳窗口确定后,引入参与建模的遥感特征变量作为自变量,实测森林蓄积作为因变量,引入多元逐步回归分析模型、随机森林模型进行建模,并根据输出的回归分析中决定系数(R2)、均方根误差(RMSE)、相对均方根误差(rRMSE)、容差等参数对模型拟合效果进行评价。

3 结果与分析

3.1 多元逐步回归模型

通过分析候选的光谱特征、纹理特征和蓄积量的相关性,从中挑选出相关性显著的特征变量进行组合,采用多元线性回归方法,建立森林蓄积量与各个光谱因子、纹理特征因子之间的多元回归模型。随着引入特征变量的增加,一方面提高了建模精度,另一方面也会引起共线性问题。当设置13×13的窗口,筛选引入特征变量为B3(红光波段),DI2(第二波段相异性),EN2(第二波段熵),SM2(第二波段二阶矩阵),CO3(第三波段相关性)时,拟合效果较好,R2增加到0.722,RMSE减少至36.305,容差为0.226,拟合效果如表3所示。

表3 不同因子多元逐步回归模型拟合效果

本文基于3×3,5×5,7×7,9×9,11×11,13×13,15×15,17×17,21×21等9种不同窗口提取的遥感特征变量与森林蓄积量进行多元逐步回归分析。实验表明:筛选引入特征变量为B3,DI2,EN2,SM2,CO3,在3×3至13×13窗口区间,容差均大于0.1,R2呈递增趋势,RMSE随之递减；在13×13至21×21窗口区间,容差均大于0.1,R2呈递减趋势,RMSE随着递增。因此,在13×13窗口下,拟合精度最优,R2(0.722) 最高,RMSE(36.305m3/hm2)最小,拟合效果如表4所示。

表4 不同窗口多元逐步回归拟合效果

遥感影像因分辨率不同对应像元面积也存在差异,GF-1单波段遥感影像单个像元面积为2m×2m,在遥感特征变量提取过程中,以13×13的移动窗口取平均值,代表676 m2遥感特征值,与样地规格半径为14.57m,面积667m2实测的圆形样地面积尺度上可较好匹配,经过上述特征因子、窗口比选,多元逐步回归森林蓄积量最优估测模型为:

y=-570.884+0.131x1-0.054x2+3.274x3+8.814x4-0.002x5

式中:x1为DI2；x2为B3；x3为EN2；x4为SM2；x5为CO3;y为样地森林蓄积量。

3.2 随机森林建模

本文提取的24个特征变量(表5),分别为光谱特征、纹理特征、归一化植被指数、比值植被指数、差值植被指数,由特征变量的重要性大小排序(图2)可知,光谱特征中B3(红光波段)、纹理特征中SM3(第三波段二阶矩阵)、DI2(第二波段相异性)、HO2(第二波段协同性)、ME3(第三波段均值) SM2(第二波段二阶矩阵)、DI3(第三波段相异性)、EN2(第二波段熵)对蓄积估测影像较大。mtry和ntree为影响随机森林建模效果2个较为重要的参数。mtry表示平均树深,一般情况,默认设置为全部自变量数量的1/3；ntree表示模型默认会生成多少株树,ArcGIS pro的默认值为100,本研究选择100。

表5 遥感特征变量相关系数分析表

图2 特征变量重要性排序

随机森林模型精度也受窗口大小的影响,对比3×3,5×5,7×7,9×9,11×11,13×13,15×15,17×17,21×21等9种不同窗口提取的遥感特征变量与森林蓄积量进行随机森林建模分析。实验表明:当窗口设置为13×13时,模型的拟合效果最优,具有最大值R2(0.955),最小值RMSE(16.305 m3/hm2)。相同尺度窗口下,对比多元逐步回归模型,随机森林算法蓄积量估测模型拟合效果更好(表6)。

表6 不同随机森林蓄积估测建模拟合效果

3.3 模型验证

实际调查马尾松样地点共180个(其中训练样本120个,验证样本60个),其中训练样本与验证样本相互独立,对比多元逐步回归、随机森林模型验证样本中实测值与预测值的散点图拟合效果,2个模型的决定系数R2分别为0.722和0.955,随机森林模型散点分布更靠近预测趋势线,分布规则性更强,因此随机森林模型拟合效果更优(图3)。

图3 不同预测模型散点分布图

由上述模型中精度检验指标得出:随机森林模型精度优于多元逐步回归模型,多元逐步回归模型参数(R2=0.722,RMSE=38.77,rRMSE=35.33%)；随机森林模型参数(R2=0.955,RMSE=17.31,rRMSE=14.15%),检验结果如表7所示。

表7 不同模型精度检验结果

4 结论与讨论

1)优化了特征变量筛选组合方法。本文将提取的遥感特征因子作为特征变量,根据随机森林建模中输出变量的重要性进行排列,选取对建模贡献较大的因子组合,采用多元逐步回归方法,筛选拟合精度较优、共线性和数据冗余问题较小的因子作为特征变量。该方法减少了多元逐步回归中因子组合迭代次数,也在一定程度上解决了因子间共线性和数据冗余问题,在因子筛选过程中发现纹理特征对森林蓄积建模精度更为敏感,最终选取DI2,B3,EN2,SM2,CO3作为建模特征变量。

2)有效解决了遥感特征提取单元与地面监测尺度匹配性问题。对比多元逐步回归方法和随机森林算法,构建3×3,5×5,7×7,9×9,11×11,13×13,15×15,17×17,21×21等9种窗口下蓄积估测模型,分析各窗口建模拟合精度指标,发现13×13窗口下拟合精度指标最优,且13×13的移动窗代表676 m2遥感特征值,与实测的圆形样地面积667m2尺度上匹配性高。因此，寻求特征变量提取单元与地面监测单位对应最优尺度是本研究的创新点。

3)随机森林模型拟合效果优于多元逐步回归模型。基于不同遥感特征变量、不同窗口构建蓄积量估测模型,相比于多元逐步回归模型,随机森林建模的拟合效果、精度检验更优,适应性更强。选取DI2,B3,EN2,SM2,CO3作为建模特征变量,以13×13窗口建立蓄积估测模型,拟合效果和检验精度分别为多元逐步回归模型(R2=0.722,RMSE=38.77,rRMSE=35.33%)、随机森林模型(R2=0.955,RMSE=17.31,rRMSE=14.15%)。因此,选取适宜窗口、较为敏感的特征变量,随机森林模型可提高建模精度。

4)本文主要考虑了光谱、纹理特征因子进行组合,而地形、环境因子等未参与建模分析。研究主要以GF-1遥感影像作为数据源,其它不同分辨率遥感影像对应的最佳建模窗口尺度需继续深入分析研究。