基于GF-1PMS影像的森林郁闭度定量估测模型研究

魏智海，张志超，李崇贵*

(1.国家林业和草原局 西北调查规划设计院，陕西 西安 710048；2.西安科技大学 测绘科学与技术学院，陕西 西安 710054)

郁闭度是森林资源调查中的重要林分因子，与此同时，它也是森林生长状况和结构的重要指标，对制定森林经营管理方案和森林总体规划设计研究具有重要意义[1-3]。传统森林郁闭度调查方法主要有目测法、样线法、样点法和树冠投影法等，但随着林业工程建设的稳步推进，传统的郁闭度调查方法不仅耗时费力，而且获取的数据量非常有限，无法满足大范围内郁闭度研究的需求[3]。

近年来，伴随着卫星技术的快速发展，利用遥感影像进行森林郁闭度定量估测已经成为一种新途径。李崇贵等[4]基于Landsat TM影像采用最小二乘估计及岭估计方法对以RS和GIS为基础的森林郁闭度估测，研究发现当影响郁闭度估测的变量间存在复共线性时，岭估计具有较高精度；陈崇等[2]基于Landsat TM影像进行亚热带森林郁闭度遥感估算，研究发现利用BP神经网络建立的森林郁闭度优于逐步回归模型；谭炳香[5]利用星载EO-1 Hyperion高光谱遥感数据进行森林郁闭度定量估测，研究发现高光谱数据精细的光谱信息非常有利于提高植被参数的估测和反演精度；柏艺[6]基于3种不同数据源(Pleiades-1、SPOT-5和Landsat-8影像)，通过建立偏最小二乘法对川西常绿阔叶林的森林郁闭度进行估测，研究发现对于川西南常绿阔叶林郁闭度的估测上，Pleiades-1影像更具有优势；吴飏等[1]采用逐步回归方法，结合光谱、地形和纹理特征进行郁闭度遥感估测，研究发现结合图像纹理特征的方法比传统的只基于光谱或地形特征的方法在估测精度上有很大提高；胡振华等[7]基于Hyperion数据采用偏最小二乘法对香格里拉森林郁闭度进行遥感估测，研究发现基于样地数据郁闭度变化敏感性分析模型优于基于Hyperion影像的分段主成分分析模型；李擎等[3]基于GF-2遥感影像采用逐步回归分析对天山云杉林郁闭度定量估测，研究发现基于纹理特征、光谱信息、地形因子相结合所建立的模型效果最好。Pu等[8]利用小波变换的Hyperion影像，基于逐步回归方法估测森林郁闭度，研究表明Hyperion影像估测精度可以达到近85%。然而，这些方法未对不同样地数据适用的模型进行分析，并且未考虑不同的森林覆盖类型对于郁闭度估测精度的影响。

因此，本研究以内蒙古自治区东北部柴河林业局为研究区域，综合利用该区域的GF-1遥感影像、二类调查小班数据、DEM数字高程模型以及林地一张图数据，采用k-NN、稳健估计以及偏最小二乘法建立以样地为单位的森林郁闭度估测模型，辅之以二类调查小班数据和现地实测数据进行模型精度检验并挑选森林郁闭度最优估测模型。

1 材料与数据

1.1 研究区概况

研究以位于内蒙古自治区东北部柴河林业局为研究对象，范围在120°36′-121°46′E， 47°05′-47°45′N(图1)。柴河林业局属寒温带大陆季风气候，年平均气温约为3 ℃，年平均降水量约为550 mm。柴河林业局地势西高东低，地貌以中低山地貌为主，西部属于中山地貌，东部属于低山丘陵地貌[9]。该地区森林覆盖率>75%，西部主要兴安落叶松覆盖为主，东部则以蒙古栎、黑桦占绝对优势[10]。

图1 研究区与样地位置图像Fig.1 Study area and locations oof sample sites

1.2 数据获取与处理

1.2.1 遥感数据来源与预处理 遥感数据为通过中国资源卫星应用中心网站下载的一景GF-1遥感影像，传感器及轨道号为：GF1_PMS1_E121.5_N47.4_20160920_L1A0001840656_PMS1。该卫星上搭载了2台2 m分辨率的全色/8 m分辨率的多光谱高分辨率相机[11]，并且它是国家高分辨率对地观测系统中的首发星。因此，研究选用GF-1遥感影像的森林郁闭度定量估测对促进国产遥感卫星在森林资源调查和监测中的应用具有重要意义。在数据预处理方面，利用ENVI5.3SP1对影像进行辐射定标、大气校正、正射校正、影像融合以及几何精校正等预处理，最后通过叠加研究区域行政界线图进行裁剪。其中大气校正使用FLAASH大气校正，正射校正使用GF-1数据自带的RPC文件和30M分辨率的DEM数据进行，几何精校正以17级谷歌卫星图(无偏移)为基准。

地形特征数据为地理空间数据云(www.gscloud.cn)下载的DEM数据，其空间分辨率为30 m，该数据主要用于提取研究区的坡度(Slope)和坡向(Aspect)，亦可作为地形特征参与郁闭度定量估测。对研究区2017年林地“一张图”数据进行重分类，并通过ArcGIS生成有林地分布图(图2)。2019年9月在该影像范围内均匀选取具有代表性的47块小班进行外业现地实测，通过实测获取到的小班郁闭度调查结果参与估测精度评价和实际应用效果讨论。

图2 研究区林地分布Fig.2 Distribution of forest land in the study area

1.2.2 地面数据来源 用于研究的地面样地数据采用2016年全国森林资源二类调查小班数据。该数据通过在研究区内按照4 km×4 km的间距布设解译样地，解译样地为0.06 hm2的长方形(10 m×60 m)样地，通过将解译样地所落小班的郁闭度拟定为样地郁闭度，最终获得有林地样地54个。

2 研究方法

本研究依据样地的基础信息(位置、大小以及形状)，通过叠加到处理好的GF-1影像提取得到样地大小的影像，通过计算遥感信息、纹理信息和DEM提取的坡度、坡向等地形信息最终获取到样地大小的自变量，通过代入3种估测模型(稳健估计、偏最小二乘法和k-NN法)计算每一个样地的林分郁闭度估测值[12-14]，以此类推，最终获得整个柴河林业局的森林郁闭度空间分布专题图。

2.1 稳健估计法

稳健估计法不仅能够有效消除或者减弱异常样地对森林郁闭度估测模型的干扰能力，而且能够有效剔除含有异常值的数据，从而使得回归模型参数具有一定的“抗干扰性”和“稳定性”[13,15]。人工有规则地挑选样地会使得模型无法满足普适性原则，因此，样地挑选应遵循随机性原则。本研究采用随机抽样的方式进行样地抽取，并且采用m：n-m比例抽取训练数据集和验证数据集。假设样地郁闭度为v，自变量个数为q，则

v=β1x1+β2x2+…+βqxq+β0

(1)

故m个样地郁闭度V与q个自变量之间的关系可以表示为

V=Xβ+eV,(eV)=0,Cov(eV)=σ2I,V=(v1,v2,…,vm)′,eV=Xβ-V

(2)

稳健估计需要满足式(3)要求才能使用，即要求误差的绝对最小值作为极值函数。

(3)

式中:ρ(ei)是极值函数，采用迭代函数法进行求解，具体步骤及原理详见文献[15]。

2.2 偏最小二乘法

偏最小二乘法是一种新型的回归方法。一方面，集成了主成分分析的优点，能够对多个自变量和因变量分别进行分析；另一方面，采用了相关性最大化等原则进行模型建立，可以有效解决变量多重相关性造成的不利影响，提高所模型模的可靠性[16-17]。具体步骤可参考文献[15,18]。

2.3 k-最近邻法

k-NN法是一种简单的机器学习算法，同时也是比较典型的非参数化方法。该方法是通过找出一个样本的k个最近“邻居”并将其属性按照取平均的方式赋值给该样本以此得到样本属性。同时，k-NN法与其他建模方法相比较，该建模方法不受自变量间不存在复共线性的限制[18-19]。公式如下

(6)

(7)

式中:Wpi,q为权重矩阵，本研究选用反距离函数计算得到；pi为训练样地点，q为待估测点；Dpi,q为待估测点q和训练样地点pi之间的距离，本研究选用欧氏距离进行计算[18-20]；t为距离分解因子，本研究取t=2；upi为点pi的森林郁闭度；upi为待估测点q的森林郁闭度。

采用k-NN法建模时，选用交叉验证的方式消除研究区样地数量少的误差，最终通过比较RMSE的方式优选最佳k值[13,20-21]。RMSE(RMSE)计算公式如下：

(8)

3 结果与分析

3.1 自变量设置与变量优选

研究发现[1,7,3]，遥感、地形和纹理特征都会影响森林郁闭度估测，因此本研究综合考虑这3个方面特征。地形特征包括坡度和坡向(阴坡和阳坡)；遥感特征包括影像波段灰度信息和比值波段。影像波段灰度信息包括蓝(B)、绿(G)、红(R)、近红外(Nir)4个波段的灰度信息；比值波段主要包括归一化植被指数(NDVI)、增强型植被指数(EVI)、差值植被指数(DVI)、修正植被指数(mNDVI)、简单比值植被指数(SR)、GF42。比值波段构建公式如表1所示。

表1 比值波段构公式Table 1 Ratio band structure formula

本研究选用的纹理特征信息主要包括均值(Mean，ME)、方差(Variance，VAR)、协同性(Homogeneity，HOM)、对比度(Contrast，CON)、相异性(Dissimilarity，DIS)、熵(Entropy，ENT)、二阶矩(Angular Second Moment ，ASM)和相关性(Correlation，COR)8种信息。现有研究表明，不同窗口大小的纹理信息对郁闭度估测的影响不同，而且某些纹理信息之间存在显著的相关性，不同纹理信息之间存在大量冗余信息[1,3]。本研究首先对GF-1影像进行降维处理，提取其第1主成分；然后利用ENVI5.3SP1提取第1主成分图像在不同窗口大小下的纹理图像，窗口大小分别设置为3×3、5×5、…、41×41，方向设置为135°，步长设置为1，灰度级压缩至64级。分别以提取的3×3窗口至41×41窗口下的8个纹理信息为自变量，以样地郁闭度为因变量，利用SPSS分别进行多元逐步回归分析，所得模型精度如图3所示。

通过模型的R、R2以及RMSE来衡量不同窗口的纹理信息对研究区森林郁闭度定量估测的贡献[22]。分析图3可知，当窗口大小为3×3时，R和R2均最大，RMSE最小，该窗口大小的纹理信息对森林郁闭度估测有较大的贡献。

图3 基于不同窗口纹理信息的逐步回归分析模型精度Fig.3 Model accuracy of stepwise regression analysis based on texture information of different windows

同一窗口下不同纹理特征之间存在显著的相关性(表2)。对3×3窗口下的各纹理特征信息以及各纹理特征信息与样地郁闭度进行相关性分析，不难发现，对比度与方差，相异性与协同性、对比度，对比度与熵，二阶矩与熵之间的相关性均>0.9，存在显著的相关性，综合考虑各纹理特征信息之间以及各纹理特征信息与郁闭度之间的相关性，最终选取均值、方差、协同性、二阶矩以及相关性共5个纹理特征作为自变量因子，在最大限度保留有效信息的基础上减少冗余信息对建模的影响。

表2 相关性分析Table 2 Correlation analysis

综上所述，本研究从遥感、地形以及纹理特征3个方面共得到15个原始特征，但过多的变量会造成数据冗余，使计算量大为增加，并且影响郁闭度估测模型的预报能力。本研究采用平均残差平方和准则和穷举算法进行建模变量优选[16,18,23]。平均残差平方和表达式为

(9)

式中:q为变量的个数，Y为因变量矩阵，X为自变量矩阵。

3.2 模型构建

为研究自变量因子设置不同对森林郁闭度估测的影响规律，本研究以遥感信息为基础，通过对4种方案进行建模变量优选，最终确定分3种方案设置可能影响郁闭度估测的自变量因子。方案1：遥感信息；方案2：遥感和地形信息；方案3：遥感、地形和纹理信息。对于每种方案，稳健估计和偏最小二乘法均经过变量优选获取最终参与建模的自变量因子，k-NN法使用方案内的所有自变量因子进行建模。每种方案均进行35次随机抽样建模，综合分析所有建模结果从而减少随机抽样给建模带来的偶然性影响。各方案建模效果如表3所示。

表3 3种方案建模综合效果比较Table 3 Comparison of the comprehensive effects of three schemes of modeling %

分析表3，横向比较可知，3种建模方案中k-NN法的模型精度标准差最小，故该方法建模效果是最好的；而相较于k-NN法而言，稳健估计和偏最小二乘法的标准差均较大并且两者建模效果相近，建模效果容易受到随机抽样的偶然性影响。纵向比较可知，在遥感信息的基础上增加地形信息或纹理信息作为自变量因子能够有效提升参数化方法的建模效果，并且标准差减小，建模稳定性增加，但k-NN法的建模效果前后没有发生明显的变化。

为研究3种建模方法在森林郁闭度定量估测中的实际应用效果，在方案3中本研究共对样地进行35次随机抽样，抽样的原则按照模型残差相对均方根误差与预报偏差相对均方根误差相近且两者都在1倍标准差范围内，而后对模型结果进行筛选，最终分别挑选出3种方法中的最具代表性的郁闭度估测模型。各模型精度如表4所示，参数化方法的最优变量为：红光波段(R)、方差(VAR)、二阶矩(ASM)、EVI、GF42、mNDVI、NDVI、SR、坡度(S)。

表4 各模型精度Table 4 Accuracy of each model (%)

由图4可见，对郁闭度估测模型中的k取不同值时发现，基于k-NN法的郁闭度估测模型的最佳k值为24。

图4 郁闭度估测模型K取不同值时均方根误差Fig.4 The root mean square error when the canopy closure estimation model K takes different values

采用稳健估计法所建模型的估测方程如下

采用偏最小二乘法所建模型的估测方程如下

3.3 估测精度分析

利用GF-1PMS影像、坡度和坡向等地形信息以及纹理信息，获取各样地大小单元对应自变量的值，建立以样地大小为单位的郁闭度估测模型，最终获得各单元对应的郁闭度[13-14]。最后利用均值法对郁闭度估测值进行统计，获取小班郁闭度。分别以47个现地实测小班数据和随机抽取的227个二类调查小班数据(约占研究区有林地小班的5%)为标准评价各郁闭度估测模型的实际应用效果。最终3种郁闭度估测模型的实际应用效果如表5所示。

表5 各郁闭度估测模型的实际应用效果Table 5 The actual application effect of each canopy closure estimation model %

依据表5可知，非参数化估测方法的实际应用效果优于参数化估测方法，整体而言3种郁闭度估测方法的应用效果均能满足森林资源二类调查的实际应用需求，但k-NN法和稳健估计的实际应用效果明显优于偏最小二乘法。

4 结论与讨论

森林郁闭度定量估测是森林遥感估测研究中的一项主要内容[24]，传统的方法未对不同样地数据适用的模型进行分析，未提取出针对不同样地数据最优的模型，估测精度难以提高。本研究对传统方法中存在的问题，运用GF-1卫星遥感影像提取遥感信息以及纹理信息，结合样地数据和地形信息，建立以样地大小为单位的森林郁闭度估测模型，最终反演获得研究区森林郁闭度，并且模型估测精度有很大提高。

稳健估计、偏最小二乘法和k-NN法3种郁闭度估测方法的应用效果均能满足实际需求，但k-NN法无论是以现地实测小班数据为标准还是以二类调查小班数据为标准，估测精度均达到85%以上，具有较好的估测效果，明显优于其他2种参数化方法。这说明GF-1 PMS 多光谱影像在森林郁闭度定量估测方面具有一定的潜力。

2种参数化方法的模型精度相近，但实际郁闭度估测效果却存在较大的差别，稳健估计的实际估测效果明显优于偏最小二乘法。

关于估测稳定性评估方面，k-NN法在森林郁闭度遥感估测中稳定性较好，而稳健估计法和偏最小二乘法估测模型不稳定，估测效果波动性较大。

在遥感信息的基础上增加地形信息或纹理信息作为自变量因子能够有效提升参数化方法的建模效果，增加建模稳定性，但对于k-NN法并没有明显的影响。遥感信息完全可以满足利用k-NN法进行森林郁闭度估测的需求，因此在进行郁闭度估测时，应根据选用的方法选择合适的建模变量因子组合。