林业科学  2016, Vol. 52 Issue (1): 18-29   PDF    
DOI: 10.11707/j.1001-7488.20160103
0

文章信息

曹宇佳, 陈尔学, 李世明
Cao Yujia, Chen Erxue, Li Shiming
基于多源数据的省级树种(组)成数空间分布信息估测方法
Estimation of Provincial Spatial Distribution Information of Forest Tree Species (Group) Composition Using Multi-Sources Data
林业科学, 2016, 52(1): 18-29
Scientia Silvae Sinicae, 2016, 52(1): 18-29.
DOI: 10.11707/j.1001-7488.20160103

文章历史

收稿日期:2015-03-04
修回日期:2015-05-25

作者相关文章

曹宇佳
陈尔学
李世明

基于多源数据的省级树种(组)成数空间分布信息估测方法
曹宇佳, 陈尔学 , 李世明    
中国林业科学研究院资源信息研究所 北京 100091
摘要[目的] 利用能够反映植被季相变化和物候差异的中空间分辨率高重访周期遥感数据以及其他多源数据,提取区域树种(组)成数空间分布信息,间接表达主要树种(组)的空间分布,为大区域树种(组)空间分布制图提供新的方法和思路。[方法] 以吉林省为试验区,以250 m空间分辨率的MODIS NDVI 8天合成时间序列数据和国家森林资源连续清查固定样地数据为主要数据源,综合利用气象观测数据和地形数据,基于梯度最近邻(GNN)方法对省级树种(组)成数进行估测。首先利用典型对应分析(CCA)对特征变量进行特征变换;然后采用k-NN方法对树种(组)成数进行分层估测,并对k-NN方法中的k值进行优选,分析k-NN估测精度随k值的变化规律;最后基于9个县的森林资源二类调查样地和省级一类清查固定样地数据,对树种(组)成数分布图进行精度检验。[结果] 对在吉林省分布较广的蒙古栎、白桦、紫椴、春榆、杨树、胡桃楸和长白落叶松7个树种(组)成数进行估测,并制作相应的树种(组)成数空间分布图。估测结果表明,树种成数分布与固定样地成数分布呈现出一致的空间分布特征。其中,县级尺度下的k-NN预测精度检验结果为:R2为0.83,RMSE为0.35;在20 km×20 km,30 km×30 km,40 km×40 km和50 km×50 km 4个尺度下的k-NN估测结果显示,各类树种(组)在40 km×40 km和50 km×50 km尺度下的估测结果较优,春榆在各个尺度下的估测精度均较高,其平均RMSE为0.35,蒙古栎的估测精度相对较低,其平均RMSE为0.65。在不同尺度下的估测结果表明,随着k值的增加,RMSE均呈现先快速减小、后趋于相对平衡的趋势,根据该规律可确定最佳k值。另外,k-NN分层估测的估测精度高于k-NN直接估测的估测精度,其在不同尺度下的RMSE相对直接估测的结果均低0.1左右。[结论] 本文提出的基于多源数据的森林树种(组)成数空间分布估测方法是一种有效的森林参数估测方法,基于该方法能够获取较高精度的树种(组)成数空间分布图。为了得到最佳的估测效果,需要对k-NN方法中的k值进行优选,该值将随试验区和数据有所不同。另外,采用分层估测的策略可以有效提高最终估测精度。
关键词多源数据    GNN    CCA    k-NN    MODIS NDVI    树种成数    制图    
Estimation of Provincial Spatial Distribution Information of Forest Tree Species (Group) Composition Using Multi-Sources Data
Cao Yujia, Chen Erxue, Li Shiming    
Research Institute of Forest Resources Information Techniques, CAF Beijing 100091
Abstract: [Objective] Remote sensing technique provides a highly effective means for extracting tree species (group) spatial distribution information. The objective of this paper is to develop a method for estimating the provincial spatial distribution information of forest tree species (group) composition using multi-sources data. Thus it could indicate the spatial distribution information of the main tree species (group) and provide a new method for extracting vegetation information in large area. [Method] The experiments were carried out over the test site of the whole Jilin Province. The time series MODIS NDVI product of 250 m pixel size and 8 days cloudy free composite and the permanent forest plot data collected by the national forest inventory (NFI) were used as the key data sources. The weather observation data and topography data were also integrated into the data sources. We developed a gradient nearest neighbor (GNN) based approach for estimating provincial forest tree species (group) composition distribution information. Firstly, the method of canonical correspondence analysis (CCA) was implemented to extract effective composited features from the original dataset. Secondly, the k-nearest neighbors (k-NN) method was applied on the extracted feature space to estimate forest tree species (group) composition number using one two-layer stratification scheme. As the value of k needs to be determined, the changing trend of k-NN estimation accuracy with the k values was analyzed. Finally, the estimation accuracy for each tree species (group) of the developed method was validated using the forest plot data of 9 counties collected by the forest resources inventory in second level and the forest plot data collected by the NFI as reference. [Result] 7 tree species (group) composition numbers including Quercus mongolica,Betula platyphylla,Tilia amurensis,Ulmus davidiana,Populus,Juglans mandshurica and Larix olgensis were extracted and the corresponding distribution maps were produced. The results showed a good consistency with the fixed plots in field. Taking county as statistic unit, the following quantitative technical targets have been achieved:the coefficient of determination (R2) was 0.83, and the RMSE was 0.34. Specifically, the accuracy has been further validated by dividing the whole coverage of Jilin Province into grids of 20 km×20 km,30 km×30 km,40 km×40 km and 50 km×50 km, taking the forest plot data collected by the NFI as reference and the grid as statistic unit. Better results could be achieved at the scale of 40 km×40 km and 50 km×50 km. The RMSE of Ulmus davidiana composition number was 0.35 and the RMSE of Quercus mongolica composition number was 0.65. The optimal k-value could be determined for the phenomenon that the RMSE firstly reduced and then tended steady with the rising k-value. In addition, the estimation accuracy of the two-layer stratification estimation method was higher than that of the direct estimation method. The results showed that:the average RMSE of estimating tree species (group) composition using two-layer stratification estimation method was 0.1 less than that using direct estimation method.[Conclusion] The proposed method for estimating the provincial spatial distribution information of forest tree species (group) composition using multi-sources data has proved to be an effective method to estimate forest parameters. Based on this method, the distribution map of forest tree species (group) composition numbers was successfully produced with high accuracy. The results indicated that the value of k needs to be optimized in order to obtain a better result, which varies depending on the experimental area and the selected data. In addition, the estimation accuracy could be improved effectively using two-layer stratification estimation method.
Key words: multi-data sources    GNN    CCA    k-NN    MODIS NDVI    tree species composition    mapping    

森林作为陆地生态系统的主体,在维持生态过程和生态平衡中发挥着重要作用(张煜星等,2007)。森林资源调查可为森林资源的科学管理、森林生态过程和机制的研究等提供重要的数据支撑。森林树种或树种组[下文用"树种(组)"表示]信息是我国森林资源样地调查的重要因子之一,主要包括树种(组)的类型、成数及其空间分布。树种(组)成数是某树种的蓄积量占林分总蓄积的比重,在森林资源调查中普遍采用十分法表示,比如,在所调查样地内某个树种的成数为6,表示该树种蓄积约占样地总蓄积的60%(55.00%~64.99%); 一个调查样地内通常会出现多个树种(组),所有树种(组)出现的成数之和等于10。在混交林中,树种组成系数最大(即蓄积量比重最大)为优势树种。由1个树种组成,或混有其他树种但材积都分别占不到10%的林分为纯林; 而由2个或者更多个树种组成,每种树木在林分内所占成数均少于10%的林分为混交林,其中针叶树种组成大于6成的林分为针叶混交林,阔叶树种组成6成以上的林分为阔叶混交林,而针叶或者阔叶组成占4~6成的林分为针阔混交林。目前,获取森林树种(组)信息主要采用地面样地调查的方法,不仅工效低、时效性差,而且很难制作出空间连续的树种(组)分布图。近年来,遥感以其宏观性、现势性、周期性强等优势,为森林树种(组)空间分布信息的提取提供了新的技术手段(曾庆伟,2010)。然而,当前利用遥感手段获取森林树种(组)空间分布信息的方法,依靠的主要是高空间分辨率或高光谱分辨率的遥感影像(曾庆伟等,2009; 陈尔学等,2007),该类方法虽然可以较为精准地获取森林的树种(组)信息,但相关遥感数据的获取较为困难,且成本较高,仅能在小区域内应用,无法满足林业大区域范围的应用需求。

在大区域森林资源遥感调查应用研究上,国内外对土地覆盖/利用类型(包括对林地类型的细分)的遥感分类方法研究较多,但对森林树种(组)信息提取方法的研究相对较少。近些年来,国外已经有学者利用决策树、插值、梯度最邻近等方法开展了树种(组)分类和信息提取研究,并取得了一定进展(Zhu et al.,1994; Xian et al., 2002; Ruefenacht et al., 2008)。其中,决策树分类方法对样本依赖性较大,当样本不足时分类精度较差(刘勇洪等, 2005Brus et al., 2012); 插值方法则具有空间局限性,达不到在宏观尺度上反映森林树种(组)空间分布的目的; 梯度最近邻(gradient nearest neighbor,GNN)J 算法通过估测不同树种(组)的森林参数信息,可间接得到不同树种(组)的空间分布信息,与决策树方法相比,对训练样本的依赖性较小,且不存在插值方法空间局限性的缺陷,已成功应用于森林树种(组)空间分布信息估测研究中(Wilson et al., 2012)。

GNN树种(组)空间分布信息估测是近几年出现的最新方法,该方法综合采用了典范对应分析(canonical correspondence analysis,CCA)和k-NN(k-nearest neighbors)2种非线性统计分析方法; 所采用数据综合利用了粗空间-高时间分辨率遥感信息(标准化植被指数,NDVI)、国家森林资源连续清查固定样地调查数据、气象台站观测数据和地形数据(由数字高程模型DEM提取); 其最终估测结果是每个树种(组)森林参数信息相对值的空间分布图。然而到目前为止,国内尚未见综合利用多源数据估测较大区域(如一个省全覆盖、全国覆盖等)森林树种(组)空间分布的研究报道。

以吉林省为试验区开展的以GNN为基本估测方法的多源数据省级树种(组)成数空间分布估测方法研究有以下几个特点:1)国外学者采用GNN方法分树种(组)估测相对胸高断面积的空间分布,而本文直接估测树种(组)成数的空间分布,有利于了解GNN方法对不同目标参数估测的适用性; 2)在对时间序列估测特征的降维处理上,国外学者采用基于傅里叶变换的方法,本文则采用经典的主成分变换方法,旨在寻求更简单、更容易理解和应用的方法; 3)在对估测结果的验证上,本文不仅采用国外学者基于不同尺度网格的验证方法,而且采用更加客观的独立样本检验方法,即基于二类调查样地数据的精度检验方法。

1 研究区概况与数据来源 1.1 研究区概况

吉林省(121°38'—131°19'E,40°52'—46°18'N)位于我国东北地区中部(图 1),面积为18.74 万km2,占全国总面积的2%。地势东南高、西北低,地形复杂,以中部大黑山为界,可分为东部山地和中西部平原两大地貌区。所处地理位置、地形条件和大气环流决定吉林省属于2种气候带:温带季风气候和温带大陆性气候,夏季高温多雨,冬季寒冷干燥。吉林省是我国的重要林业基地,森林主要分布在东部地区,占林地面积的94%。该省主要树种有长白落叶松(Larix olgensis)、蒙古栎(Quercus mongolica)、白桦(Betula platyphylla)、杨(Populus)、紫椴(Tilia amurensis)、胡桃楸(Juglans m and shurica)等,主要森林类型为阔叶混交林、针阔混交林和栎类林(王新闯等,2011)。

图1 研究区(吉林省)地理位置(左图)及所采用的一类清查样地空间分布(右上)和DEM(右下) Fig.1 Location of the test site ( Jilin Province) ( left) and the spatial distribution map of NFI permanent forest plots ( right up) and the DEM ( right down)
1.2 数据来源 1.2.1 样地数据

1)国家森林资源连续清查固定样地数据 吉林省第7次森林资源连续清查(简称一类清查)固定样地数据,调查时间为2005年。样地为方形,大小为0.06 hm2。本研究使用了每个样地的地理坐标、优势树种和各树种成数组成。样地的树种组成采用10成法表示,如样地内70%为红松(Pinus koraiensis)、20%为白桦、10%为春榆(Ulmus davidiana),则该样地的树种成数组成为"7红松2白桦1春榆";由于红松成数占优势,所以该林分的优势树种为"红松"。

由样地数据统计发现,全省共有18个树种(组),包括紫椴、蒙古栎、白桦、春榆、杨树、色木槭(Acer mono)、胡桃楸、黄菠萝(Phellodendron amurense)、枫桦(Betula costata)、长白落叶松、红松、樟子松( Pinus sylvestris var. mongolica )、黑松(Pinus thunbergii)、云杉(Picea asperata)、柳(Salix)、水曲柳(Fraxinus m and shurica)和臭松(Abies nephrolepis)。某树种在某样地中成数不为0的样地称为此树种的有效样地,本研究仅对在试验区分布较广的主要树种(组)进行树种成数估测,选择有效样地占样地总数(N=3 594)20%以上的树种,共有7个,分别是蒙古栎、白桦、紫椴、春榆、杨树、胡桃楸和长白落叶松,筛选后样地总数为2 964。

2)森林资源规划设计调查固定样地数据 森林资源规划设计调查(简称二类调查)主要采用小班调查和固定样地相结合的方法,县(市)、局级二类调查固定样地在省级固定样地基础上加密布设得到。本研究采用吉林省9个林业局(图 2)共8 393个森林资源二类调查固定样地数据作为验证数据,实际用到的样地信息包括样地点地理坐标、样地的优势树种类型和各树种成数组成,树种成数组成的调查精度达到了90%以上。

图2 吉林省森林资源二类调查样地点空间分布 Fig.2 The spatial distribution map of forest management inventory plots in Jilin Province
1.2.2 时间序列NDVI数据

本研究采用美国NASA EOS数据中心提供的MODIS陆地产品:空间分辨率为250 m的8天合成表面反射率数据(MOD09Q1)。

利用MRT投影转换工具,将吉林省的数据进行拼接、裁剪、转投影和格式转换,形成8天合成数据集,再将生成的数据进行标准化植被指数(normalized difference vegetation index,NDVI)反演,并叠加生成8天合成周期的MODIS NDVI时间序列数据。全年共46个时相,分辨率为250 m。

所使用的MODIS产品数据为2005年全年数据集。在MODIS表面反射率产品波段中,波段1(620~670 nm)是红光(red)波段,波段2(841~876 nm)、波段5(1 230~1 250 nm)属于近红外(nir)波段。波段2是计算NDVI的常规波段范围,本研究选择波段1和2的反射率(ρ)生成NDVI 数据:

$NDVI=\frac{{{\rho }_{nir}}-{{\rho }_{red}}}{{{\rho }_{nir}}+{{\rho }_{red}}}$ (1)

1.2.3 时间序列气象数据

为更加准确地估测树种成数,引入时间序列气象数据。本研究选用中国气象科学数据共享服务站提供的2005年吉林省24个气象站点的每月平均气温数据(℃)和每月总降水数据(mm)。每类数据全年年均有12个时相。

1.2.4 地形及地理位置数据

本研究采用空间分辨率为90 m的SRTM DEM数据产品提取吉林省的高程(H)、复合地形指数(compound topographic index,CTI)、坡度坡向指数(slope aspect index,SAI)、地理坐标($\varphi $,φ)。

CTI是表征特定景观中水分和沉积物运移的综合地形变量,其计算公式(McKenzie et al., 1999)为:

$CTI=\ln \left(\frac{{{A}_{C}}}{\tan \beta } \right)$ (2)
式中:AC为垂直于水流方向的特定汇流面积,可通过ArcGIS水文分析模型计算; β为坡度。

SAI计算公式(Nielsen et al.,1998)为:

$SAI=\sin(aspect+225)\times \left(\frac{slope}{45} \right)$ (3)
式中:aspect为坡向(°); slope为坡度(°); SAI的取值范围为(-1,1)。

另外,本文还选用了研究区每一个像元的地理坐标($\varphi $,φ)作为自变量。采用地理位置作为自变量,主要是为了控制k-NN估测算法优先选择哪些与待估测像元在地理空间上距离最近的训练样地,避免出现所选择的近邻与待估测像元处于差别很大的生态区。

提取得到90 m的H,CTI,SAI,$\varphi $,φ栅格数据后,再采用三次卷积法将其重采样为250 m像元大小。

2 树种(组)成数估测方法 2.1 时间序列特征变量处理 2.1.1 MODIS NDVI时间序列数据滤波处理

由于MODIS合成数据采集过程中会受到太阳高度角、观测角度以及云、水汽和气溶胶等因素的影响,使得时间序列数据出现异常值或者缺失数据(穆少杰等,2012),因此,为了使NDVI时间序列数据能正确反映植被真实的季节性变化规律,需要对数据进行滤波处理。本研究采用Savtzky-Golay(S-G)滤波法对NDVI时间数据进行平滑滤波处理。

2.1.2 MODIS NDVI时间序列数据降维处理

MODIS NDVI时间序列数据全年46个时相波段之间存在高相关性和高冗余度,因此本研究通过主成分分析(principal component analysis,PCA)技术对滤波后的MODIS NDVI时间序列影像进行降维处理。PCA可以将一组彼此相关的变量变换为一组新的相互独立且相互正交的变量,有效降低高维度数据的冗余信息,减少后续分析的计算量和复杂度,同时又尽可能地避免信息的丢失。本研究选用PCA降维处理后的前5个主成分分量作为自变量用于树种(组)成数的估计。

2.1.3 月平均气温和月总降水量数据的插值与降维处理

吉林省气象样地点随机分布,且气象值呈正态分布,符合克里金插值要求,因此,本研究采用克里金插值法将气象样点数据插值成空间分辨率为250 m的栅格数据。最终分别得到全年12个月的月平均气温和月总降水栅格数据。

为去除气温与降水时间序列数据中的冗余信息,运用PCA方法对全年月平均气温和月总降水量数据进行降维处理。最终选用气温和降水的前3个主成分分量作为自变量用于树种(组)成数的估计。H B ┣┣三级标题┫┫B 2.1.4 树种(组)成数估测因变量和自变量 H 将所有自变量数据重投影为Albers Conical Equal Area,WGS84,空间分辨率重采样为250 m。

设响应变量(Y)为样地或一个250 m×250 m像元内优势树种(组)的成数。在估测模型建立阶段,Y就是一类清查样地7个树种的成数。由于是分别对树种(组)建立估测模型,因此每个样地的响应变量就是该样地某树种的成数(若该树种未在该样地中出现,Y就设为0),是一个一维变量。

自变量(X)由PCA变换的NDVI时间序列数据、PCA变换的月平均气温数据和月累计降水数据、地形数据(包括DEM、CTI和SAI和森林植被的空间位置)组成。

2.2 GNN估测方法

GNN方法是将直接梯度排序分析(如CCA)与最邻近算法(k-NN)相结合的一种算法(Ohmann et al., 2002),主要分为4个步骤:1)通过CCA分析样地Y变量与X变量之间的定量关系,得到新的特征变量组; 2)根据每个像元中特征变量组的特征值,选择最佳特征变量组合; 3)根据最邻近算法原理确定待测像元周围的k个邻近样地点; 4)以k个邻近样地点的Y值估算待测像元的Y值,样地与待测像元之间的权重采用欧氏距离计算得到。

2.2.1 CCA

CCA是基于对应分析发展而来的一种排序方法,将对应分析与多元回归分析相结合,每步计算均与自变量进行回归,又称多元直接梯度分析(张元明等,2004)。其基本思路是,在对应分析的迭代过程中,每次得到的Y变量均与自变量进行多元线性回归,即:

${{Y}_{j}}={{b}_{0}}+\sum\limits_{k=1}^{q}{{{b}_{k}}}{{U}_{kj}}$ (4)
式中:Yj为第j个样点的值; b0为截距(常数);bk(k=1,2,…,qq为自变量个数)为样点与第k个自变量之间的回归系数:Ukj为第k个自变量在第j个样点中的测量值。

2.2.2 k-NN

k-NN是一种典型的非参数方法,基于观测点和预测点之间的空间相似性关系进行单变量或多变量预测(Franco-Lopez et al.,2001; McRoberts et al.,2002; Tomppo,1991)。利用k-NN算法进行参数估计的基本原理为:记p为待测点,pi为邻近参考点,Dpi,p为两点之间的光谱距离,其中,参考点pi的树种成数是已知的。Dpi,p用于衡量样本之间的相似度,其距离值越小,表明相似度越大,反之则表明相似度越小。对于待测点p,找出其光谱空间最近邻的k个样地点p1p2,…,pk,其中,Dp1,p≤Dp2,p≤…≤Dpk,p。待测像元p的树种成数(Wp)通过k个参考点相应的树种成数加权平均获得,如式(5):

${{W}_{p}}=\sum\limits_{i=1}^{k}{{{\omega }_{{{p}_{i}},p}}}\times {{W}_{{{p}_{i}}}}$ (5)
其中,权重值通过计算距离的倒数获得,如式(6):
${{\omega }_{{{p}_{i}},p}}=\frac{1}{{{D}_{{{p}_{i}}}}}/\sum\nolimits_{j=1}^{k}{\frac{1}{D_{{{p}_{j}},p}^{t}}}$ (6)

k-NN实质上是一个常用于空间插值的反距离加权平均法,当k=1时,k-NN即为最邻近距离法。为保证k-NN算法估测树种成数的总体精度,需要对k值进行优选。k为距离分解因子(distance decomposition factor),k越大,估计结果越容易受光谱距离近的值影响,k一般取0,1,2。参考点和目标点之间的光谱距离可以采用多种距离来度量,最常用的为欧氏距离(Euclidean distance)和马氏距离(Mahalanobis distance)。研究表明,在森林相似的条件下,k-NN用于森林定量估测时,欧氏距离和马氏距离相差不大(McRoberts et al., 2002)。本研究选用欧氏距离,如式(7):

${{D}_{{{p}_{i}},p}}=\sqrt{\sum\limits_{i=1}^{n}{({{X}_{{{p}_{{}}}}}-{{X}_{p}})}}$ (7)
式中:xp表示待测像元p的光谱向量;xpi表示pi所在像元的光谱向量。

2.3 分层估测树种成数 2.3.1 样地树种成数与自变量的定量关系分析

基于R语言的Vegan程序包,利用CCA模型分析树种成数(Y)与自变量X之间的定量关系。首先,运用蒙特卡罗置换检验环境变量与树种成数分布相关的显著性,结果为0.001,说明排序结果可以接受环境因子对树种分布的解释量;然后,在CCA模型中将X变量与Y变量进行回归分析。最终得到7组X变量的线性组合表示与Y变量的回归关系,其中,前4组特征变量特征值的总和占总特征变量特征值的93.2%,因此本研究选用前4组X变量的线性组合作为新的特征变量。

2.3.2 k-NN分层估测树种成数

MODIS影像的像元大小为6.25 hm2,森林资源一类清查样地大小为0.06 hm2,二者相差很大; 且MODIS像元多为混合像元,样地只能表示混合像元中的某一类地物。因此需要将一类清查样地进行分类,如分为高植被覆盖度样地和低植被覆盖度样地,分层估测以减小混合像元对估测精度的影响。

相比250 m分辨率的MODIS数据,一类清查样地数据可以代表一个TM像元的真实情况。将TM影像重采样为250 m,一个MODIS影像的像元相当于TM影像的10×10个像元。基于TM提取植被覆盖度,对于森林而言,此植被覆盖度为冠层覆盖度。有研究指出,当植被覆盖度在0.35以下时,多为低植被覆盖度区域(Liu et al., 2013; Buyantuyev et al., 2007; Small et al.,2006),因此选用TM植被覆盖度值为0.35作为低植被覆盖层和高植被覆盖层的分界值。

若一类清查样地位置落在TM低植被覆盖层,则这些样地数据将作为低植被覆盖度一类清查样地; 反之,为高植被覆盖度一类清查样地。分别对2类样地进行k-NN估测,得到2个估测结果(Y1Y2)。设在10×10个TM影像像元中(等效于一个MODIS像元大小)低植被覆盖度像元的百分比为R1,高植被覆盖度像元的百分比为R2,则最终的估测值(Y)为二者的加权和:Y=R1×Y1 +R2×Y2

2.4 估测精度评价方法

由于本研究是基于一个省级单位进行森林树种(组)成数估测,估测结果主要是一个省级区域各树种(组)丰富度的大致分布情况,并不是对小区域(落实到山头地块)的树种(组)进行详细研究,因此,在进行精度检验时,不必在像元尺度下进行检验,只要在县级尺度下满足估测精度即可。

2.4.1 县级尺度下的精度验证

本研究选取吉林省9个县的二类调查样地对估测结果进行检验。将每个县级区域内所有二类调查样地树种成数的平均值和k-NN遥感预测值的平均值分别作为该区域内树种成数的实测值和预测值,计算二者的决定系数(coefficient of determination,R2)和均方根误差(RMSE)以及每一个树种的R2和RMSE:

${{R}^{2}}=1-\left[ \frac{\sum\limits_{i=1}^{n}{{{({{y}_{i}}-{{{\hat{y}}}_{i}})}^{2}}}}{\sum\limits_{i=1}^{n}{({{y}_{i}}-{{{\hat{y}}}_{i}})}} \right];$ (8)
$RMSE=\sqrt{\frac{1}{n}{{\sum\limits_{i=1}^{n}{({{y}_{i}}-{{{\hat{y}}}_{i}})}}^{2}}}$ (9)
式中:yi为实测值; ${{\hat{y}}_{i}}$为预估值; A 为实测值的平均数; n为样本数。

2.4.2 多尺度下的精度检验

由于采用二类调查样地进行精度检验,只是针对吉林省的9个县进行抽样检验,无法实现全省覆盖的精度检验,因此本文将全省划分成一定大小的网格,以落入网格内的一类清查固定调查样地数据为参考,以网格为统计单元进行精度检验。

为保证每个估测尺度统计单元内有足够多的样地用于估计响应网格的树种成数,将估测尺度划分为20 km×20 km,30 km×30 km,40 km×40 km和50 km×50 km共4个估测尺度,4个估测尺度分别对应576,271,161和114个尺度单元(图 3),图中每个方格代表一个尺度单元,选取的尺度单元至少包含5个样地。以一定尺度统计单位内多个样地树种成数的均值作为该统计单元树种成数的实测值,以k-NN遥感预测值的均值作为该统计单位树种成数的预测值,计算二者的均方根误差(RMSE)。

图3 4 个精度评价尺度下的样地分布 Fig.3 Sample plot distribution in four accuracy assessment scales
3 结果与分析 3.1 k值的优选

以白桦为例,分析白桦在4个尺度下估测精度随k值的变化(图 4)。随着k值的增加,RMSE呈现先快速减小、后趋于相对平衡的趋势。当k=6时,白桦树种成数的RMSE在4个尺度下几乎都是最低,因此选取k=6为最佳。从图 4还可以看出,在不同估测尺度下,相同k值的估测精度虽然不同,但是整体变化趋势一致。其他6个树种均按照此方法进行k值的优选,蒙古栎、紫椴、春榆、杨、胡桃楸和长白落叶松的k值选优结果分别为5,6,6,6,5和7。

图4 白桦的k-NN 估测树种成数的RMSE Fig.4 k-NN estimation RMSE of White Birth tree species composition in multi-scales
3.2 k-NN分层估测对估测精度的影响

以白桦为例,对比分析k-NN分层估测与k-NN直接估测在4个尺度下的RMSE。从图 5可以看出,k-NN分层估测的估测精度在4个尺度下均小于k-NN直接估测的RMSE,说明采用k-NN分层估测方法可提高树种成数的估测精度,减小估测误差。

图5 k-NN 分层估测与k-NN 直接估测的RMSE Fig.5 RMSE of direct and stratification estimation in multi-scale using k-NN
3.3 树种成数的估测结果

对每个树种选用估测精度最高时所对应的k值进行k-NN分层估测,得到吉林省7个树种的树种成数分布图(图 6~12)。图 6为白桦的森林资源一类清查样地点树种成数的等级空间分布图和k-NN估测树种成数等级分布图。由图 6可知,k-NN估测树种成数分布与固定样地成数分布呈现出一致的空间分布特征,白桦多分布在1~3成。对7个树种总体而言,k-NN估测的树种成数也多在1~3成,这主要是由于吉林省森林类型多为混交林。其中,7个树种中以纯林(树种成数为7成以上)成片分布的只有蒙古栎,且多分布在吉林省东北部(图 11)。吉林省森林绝大部分分布在东部,西部只有少量杨树分布,树种成数集中在4~6成(图 12)。通过k-NN估测7个树种的树种成数分布可知每一类树种的地理分布和生长丰富度情况。在南部,除白桦和杨外,其他树种均有分布。胡桃楸主要集中分布在吉林省中部,其他树种则主要分布在吉林省东部。

图6 白桦一类清查样地点的树种成数等级分布和k-NN 估测的白桦树种成数等级分布 Fig.6 The White Birch grade distribution map of NFI permanent forest plot data and tree species composition grade distribution map of White Birch in k-NN estimation method

图7 k-NN 估测的紫椴树种成数等级分布 Fig.7 Tree species composition grade distribution map of Linden in k-NN estimation method

图8 k-NN 估测的胡桃楸树种成数等级分布 Fig.8 Tree species composition grade distribution map of Manchurian Walnut in k-NN estimation method

图9 k-NN 估测的春榆树种成数等级分布 Fig.9 Tree species composition grade distribution map of Elm in k-NN estimation method

图10 k-NN 估测的长白落叶松树种成数等级分布 Fig.10 Tree species composition grade distribution map of Dahurian Larch in k-NN estimation method

图11 k-NN 估测的蒙古栎树种成数等级分布和局部放大 Fig.11 Tree species composition grade distribution map of Mongolian Oka Tree in k-NN estimation method and partial enlarged detail

图12 k-NN 估测的杨树树种成数等级分布和局部放大 Fig.12 Tree species composition grade distribution map of Aspen in k-NN estimation method and partial enlarged detail
3.4 精度检验 3.4.1 县级尺度下的精度检验

图 13为利用吉林省9个县二类调查样地得到的县级区域内k-NN遥感估测值与样地实测值之间关系的分析结果,横坐标为样地实测值,纵坐标为k-NN估测值,虚线表示理想回归线(y=x),实线表示估测值拟合的回归线。从图 13中可以看出二者几乎重合,说明在县级尺度下k-NN模型预测效果较好。k-NN预测精度检验结果为:相关系数R2达到0.83,RMSE为0.35。表 1为7个树种的R2和RMSE,从表 1可以看出,蒙古栎的R2最高,为0.92,而春榆的R2较小。

图13 树种成数实际值和预测值散点图 Fig.13 Scatter plots of observed tree species ( group) composition number with the predicted values

表1 k-NN 估测的7 个树种的树种成数在县级尺度下的R2 和RMSE Tab.1 R2 and RMSE of k-NN estimation of tree species composition for seven tree species in county scale
3.4.2 多尺度下的精度检验

表 2 为7 个树种在20 km × 20 km,30 km × 30 km,40 km × 40 km和50 km × 50 km 4 个尺度下的估测情况。由于7 个树种中每个树种的一类清查样地点空间分布不同,且树种的生态特性不同,使得不同树种对尺度表现出的敏感性也不同。从表 2 可以看出,紫椴在40 km× 40 km 尺度下RMSE 最低,为0. 33,杨树则在20km × 20 km 尺度下RMSE 最高,为0. 53;而且,随着尺度的增大,树种成数的平均RMSE 逐渐降低。在4 个评价尺度下的平均RMSE,春榆最低,为0. 35;蒙古栎最高,为0. 65。树种的估测精度主要受样地点个数和样地点分布的影响,如果样地点多且分布集中,则估测单元内的样地点多,估测误差小。

表2 7 个树种(组)在多尺度下的估测精度 Tab.2 k-NN estimation accuracy of tree species (group) composition number for seven tree species (group) in multi-scales
4 结论与讨论

本文研发了一种基于多源数据的森林树种(组)成数空间分布估测方法,其特点是综合利用时间序列遥感数据、国家森林资源固定样地调查数据、气象观测数据和地形数据,采用先进的GNN方法进行估测。研究发现:k-NN是一种典型的非参数化估测方法,为了得到最佳的估测效果必须对k值进行优选,而且最佳k值会随不同的试验区和数据有所不同; 多尺度分层估测的RMSE均小于不分层直接估测的RMSE,平均降低了0.1,说明本文提出的分层估测方法是有效的; 采用森林资源二类调查样地数据在县级尺度下对估测结果进行了精度检验,预测结果与参考数据一致性较高,R2为0.83,RMSE为0.35; 进一步将全省划分为不同尺度的网格,以一类清查固定样地数据为参考,以网格为统计单元在不同尺度下分树种对估测结果进行精度检验,各尺度平均RMSE随尺度的增大而减小,从另一个方面说明了本文方法的有效性。

本研究为估测大区域范围树种(组)的森林参数提供了一种有效方法,该方法将粗分辨率影像与森林资源一类清查样地相结合,并采用k-NN分层方法进行估测。然而,在对样本进行分层处理时,只利用了植被覆盖度这一标准,没有考虑其他因素影响,如DEM、坡度、坡向和气候等;而且,本研究利用高时间分辨率影像可获取植被物候信息的特点进行树种成数的估测,但影像空间分辨率较低。因此,如何将高时间分辨率影像与较高空间分辨率影像结合,进行大区域树种(组)成数的估测,为下一步的研究方向。

本研究在国内首次制作了省级全覆盖、空间分辨率为250 m的主要树种(组)成数分布图,并对其进行了精度检验,对国家森林资源宏观管理政策的制定、森林对气候变化响应机制的研究等具有重要的参考价值。

参考文献(References)
[1] 陈尔学,李增元,谭炳香,等.2007.高光谱数据森林类型统计模式识别方法比较评价.林业科学,43(1):84-89.
(Chen E X, Li Z Y, Tan B X, et al. 2007.Validation of statistic based forest types classification methods using hyperspectral data.Scientia Silvae Sinicae,43(1):84-89.)(1)
[2] 刘勇洪,牛铮,王长耀.2006. 基于MODIS数据的决策树分类方法研究与应用.遥感学报,9(4):405-412.
(Liu Y H, Niu Z, Wang C Y. 2006. Research and application of the decision tree classification using MODIS data. Journal of Remote Sensing, 9(4):405-412.[in Chinese])(1)
[3] 穆少杰,李建龙,陈奕兆,等.2012. 2001-2010年内蒙古植被覆盖度时空变化特征. 地理学报, 67(9), 1255-1268.
(Mu S J, Li J L, Chen Y Z, et al. 2012. Spatial differences of variations of vegetation coverage in Inner Mongolia during 2001-2010. Acta Geographica Sinica, 67(9):1255-1268.[in Chinese])(1)
[4] 王新闯,齐光,于大炮,等.2011.吉林省森林生态系统的碳储量、碳密度及其分布.应用生态学报, 22(8):2013-2020.
(Wang X C, Qi G, Yu D P, et al. 2011.Carbon storage,density,and distribution in forest ecosystems in Jilin Province of northeast China. Chinese Journal of Applied Ecology, 22(8):2013-2020.[in Chinese])(1)
[5] 曾庆伟,武红敢.2009.基于高光谱遥感技术的森林树种识别进展.林业资源管理,10(5):109-114.
(Zeng Q W, Wu H G.2009.Development of hyperspectral remote sensing application in forest species identification.Forest Resources Management,10(5):109-114.)(1)
[6] 曾庆伟.2010.基于Hyperion高光谱数据的森林类型精细识别研究.北京:中国林业科学研究院博士学位论文.
(Zeng Q W. 2010. Forest type precise identification based on Hyperion data. Beijing:PhD thesis of Chinese Academy of Forestry.[in Chinese])(1)
[7] 张煜星,王祝雄,武红敢,等.2007.遥感技术在森林资源清查中应用研究. 北京:中国林业出版社.
(Zhang Y X, Wang Z X, Wu H G, et al. 2007. Application of remote sensing technology in forest resource inventory. Beijing:China Forestry Publishing House.[in Chinese])(1)
[8] 张元明,陈亚宁,张小雷.2004.塔里木河下游植物群落分布格局及其环境解释.地理学报, 59(6):903-910.
(Zhang Y M, Chen Y N, Zhang X L. 2004. Plant communities and their interrelations with environmental factors in the lower reaches of Tarim River. Acta Geographica Sinica, 59(6):903-910.[in Chinese])(1)
[9] Brus D J, Hengeveld G M, Walvoort D J J, et al. 2012. Statistical mapping of tree species over Europe. European Journal of Forest Research,131(1):145-157.(1)
[10] Buyantuyev A, Wu J, Gries C. 2007.Estimating vegetation cover in an urban environment based on Landsat ETM+ imagery:a case study in Phoenix, USA. International Journal of Remote Sensing, 28(2):269-291.(1)
[11] Franco-Lopez H, Ek A R, Bauer M E.2001. Estimation and mapping of forest stand density, volume, and cover type using the k-nearest neighbors method. Remote Sensing of Environment,77(3):251-274.(1)
[12] Liu T, Yang X. 2013.Mapping vegetation in an urban area with stratified classification and multiple end member spectral mixture analysis. Remote Sensing of Environment, 133:251-264.(1)
[13] McKenzie N J, Ryan P J. 1999.Spatial prediction of soil properties using environmental correlation. Geoderma, 89(1):67-94.(1)
[14] McRoberts R E, Nelson M D, Wendt D G.2002. Stratified estimation of forest area using satellite imagery, inventory data, and the k-Nearest Neighbors technique. Remote Sensing of Environment, 82(2):457-468.(2)
[15] Nielsen S E,Haney A. 1998. Gradient responses for understory species in a bracken-grassland and northern dry forest ecosystem of northeast Wisconsin. Trans Wisc Acad Sci Arts Lett, 86:149-166.(1)
[16] Ohmann J L, Gregory M J. 2002.Predictive mapping of forest composition and structure with direct gradient analysis and nearest-neighbor imputation in coastal Oregon, USA. Canadian Journal of Forest Research, 32(4):725-741.(1)
[17] Ruefenacht B, Finco M V, Nelson M D, et al. 2008.Conterminous US and Alaska forest type mapping using forest inventory and analysis data. Photogrammetric Engineering & Remote Sensing,74(11):1379-1388.(1)
[18] Small C, Lu J W T.2006. Estimation and vicarious validation of urban vegetation abundance by spectral mixture analysis. Remote Sensing of Environment, 100(4):441-456.(1)
[19] Tomppo E. 1991. Satellite image-based national forest inventory of Finland. International Archives of Photogrammetry and Remote Sensing, 28:419-424.(1)
[20] Wilson B T, Lister A J, Riemann R I.2012. A nearest-neighbor imputation approach to mapping tree species over large areas using forest inventory plots and moderate resolution raster data. Forest Ecology and Management, 271:182-198.(1)
[21] Xian G, Zhu Z, Hoppus M, et al.2002. Application of decision-tree techniques to forest group and basal area mapping using satellite imagery and forest inventory data.Pecora, 15:10-15.(1)
[22] Zhu Z, Evans D L.1994. US forest types and predicted percent forest cover from AVHRR data. Photogrammetric Engineering and Remote Sensing, 60(5):525-531.(1)