近些年,我国华北雾霾天数逐渐增多,空气污染加剧. 为净化区域空气、调节气候、提高首都生态承载能力和建设绿色大都市,北京市于2012—2015年开展了百万亩平原人工造林工程. 造林完成后,亟需相关的监测方法对林木成活率、林木长势和健康、林木管护情况及造林成果进行评价. 随着对地观测技术的日渐成熟,利用遥感探查森林种类及分布情况成为经营、管理森林资源的基本手段 ^{[ 1- 2]} .

由于低、中空间分辨遥感影像在空间、纹理上的局限性 ^{[ 3]} ,使得高空间分辨率影像因结构清晰、纹理突出等优势在树种提取方面得到了广泛地应用 ^{[ 4]} . 然而,利用传统的分类算法对高空间分辨率遥感影像进行分类存在困难,Baatz等 ^{[ 5]} 针对高分辨率遥感影像的特点,考虑了地物的光谱、纹理及形状等特征,提出面向对象的遥感分类方法,该方法有效避免了基于像元分类的“椒盐效应”,极大提高了分类精度.

影像分割是面向对象分类的一个关键步骤,其结果直接关系到后续信息提取的精度,传统影像分割最佳尺度的选择通常采用试错法,通过尺度的“增长”目视分割效果以确定最优分割尺度,但遥感影像数据量大,分割效率较低且不利于推广. 许多学者提出了改进的方法,如李秦等 ^{[ 6]} 通过改进局部方差法确定了最优分割尺度;Drgut等 ^{[ 7]} 借鉴Woodcock等 ^{[ 8]} 所提出的影像局部方差(local variance,LV)和最优分辨率之间的关系,又引入了方差变化率(rate of change,ROC)来描述层次间的变化,快速计算出影像的最优分割尺度,该方法简单、实用,成为尺度选择中较为常用的方法之一 ^{[ 9]} . 分类规则的构建是面向对象分类过程中的一个热点问题 ^{[ 10]} ,研究人员一般根据区域及遥感影像的特点,人为分析和选择特征,具有一定的主观性和盲目性,选取的特征难以保证分类精度与效率. 尽管目前成熟的特征选择算法很多,但针对面向对象分类中小样本和高维特征的算法较少 ^{[ 11]} ,甄选一种快速有效的特征选择方法显得十分必要.

由于多时相遥感数据可反映不同树种的物候期,已有学者成功将其运用到树种分类. Hill等 ^{[ 12]} 指出,对于树种分类最好的时期是采用春季或秋季并加入夏季植被生长盛期影像;Tigges等 ^{[ 13]} 利用5景RapidEye影像,结合SVM分类方法,提取了柏林市城区8类树种,其结果表明,多时相影像很大程度上提高了城市树种分类的精度. 目前多时相影像主要应用于自然森林分类、城市林地植被信息或单种林地类型提取,对于树种丰富、分布复杂、类型层次多样的森林景观是否适合,则有待于进一步验证.

针对人工造林中树种混杂、景观差别大、郁闭度低等引起遥感难以高精度自动识别树种的问题,本研究利用多时相、高空间分辨率卫星影像,采用面向对象的最邻近分类法、支持向量机法及最大似然法提取区域树种分布范围,以期为平原造林成果检验、森林经营决策与规划、森林生态服务功能评价等提供基础数据和科学依据.

1 研究数据

研究区域位于北京市延庆县蔡家河流域平原造林区,该流域处于妫水河上游,是官厅水库重要水源地. 平原造林工程中,该区域先后完成人工造林面积约2300hm ²,是北京市平原地区森林面积最大的地区之一,2012年、2013年和2014年分别种植主要包括杨树( Populus tremula)、柳树( Salix matsudana Koidz)、槐树( Sophora japonica Linn)、马尾松( Pinus massoniana Lamb)、白蜡( Fraxinus chinensis Roxb)、榆树( Ulmus pumila)在内的树种. 新建林地郁闭度低,为0.1~0.5,树种多样,自然景观变化大.

选择研究范围两景Pleiades影像数据(全色波段空间分辨率0.5m,多光谱波段空间分辨率2m),时间分别为2014年5月18日、6月27日,同时获取到两景高分(GF)一号数据(全色波段空间分辨率2m,多光谱波段空间分辨率8m),时间分别为2014年3月17日、4月23日. 根据影像提供的定标系数、光谱波长信息分别对数据进行辐射校正、大气校正(Flassh模型),采用基于亮度调节的平滑滤波方法 ^{[ 14]} 对多光谱和全色影像进行融合,利用航片数据对影像进行精校正,完成数据预处理工作. 以植被长势最为茂盛的6月27日影像为基础分析数据,截取地物类型丰富2000×2000像元大小范围进行研究,其中面向对象分割和最近邻分类(K-nearest neighbor,KNN)操作借助eCognition 8.7完成,支持向量机(support vector machine,SVM)和最大似然分类(maximum likelihood classification,MLC)方法在ENVI5.0中实现.

野外调查于2014年7月—11月期间开展,由于区域为人工造林区,林木长势基本一致,选取能代表该树种典型林木2~3棵,量测其冠幅、胸径、树高、株行距,并记录GPS等信息,为后续样本的选取及精度评价奠定基础.

2 研究方法

2.1 不同地物影像特征分析

研究区域主要地物包括房屋、道路、水体、阴影、针叶树种(以马尾松为主)、杨树、榆树(金叶榆、普通榆树)、槐树、柳树、白蜡、灌木及其他树种等. 在影像上选取每类地物的典型样本不少于20个,统计其像元平均值绘制多时相影像光谱曲线图,如 图1所示,其中横坐标B、G、R、NIR代表影像蓝、绿、红、近红外波段,纵坐标值表征经大气校正后影像的像元值(反射率×10000). 从 图1可以看出4个时相中房屋、道路的光谱曲线均高于其他地物,水体光谱曲线最低. 由于光照角度及植被茂密程度的差异,4个时相中所产生的阴影大小和位置不同,导致阴影呈现出的特征趋势不一致. 从光谱曲线特征分析,夏季影像(6月27日)在4个时相中所有地物的区分效果最好,特别是金叶榆在绿波段和近红外波段上像元值呈现较高峰值,而针叶树种在初春影像(3月17日)及夏季影像近红外波段上与其他地物区分度较好,其他树种则光谱值差别较小不易区分.

为进一步分析影像、突出植被特征,采用标准假彩色方式对影像进行合成显示,4个时相影像纹理特征如 图2所示. 3月17日高分影像中以针叶树种的纹理最为清晰,呈现暗褐红色;4月23日高分影像中杨树的纹理特征较为明显,柳树等其他树种外围轮廓逐渐清晰;而5月18日Pleiades影像进一步明确了树种间的界限,且受林下地被影响较小;6月27日Pleiades影像中地被处于生长盛期,各类树种长势繁茂,金叶榆明黄特征更为明显,槐树絮状纹理突出. 因此对于仅用光谱信息不易分离的树种,多时相中纹理信息或是提高其区分度的一个重要因素.

图  1  多时相影像地物光谱曲线

Figure  1.  Mean spectral curves of land feature in multi-temporal images

下载: 全尺寸图片 幻灯片

2.2 面向对象的多尺度影像分割

多尺度分割算法实现的主要原则是所分割对象的平均异质度最小. 首先计算像元(对象)与相邻对象合并后的异质性,如果相邻对象合并后的异质性小于给定的阈值,则将其合并,对此过程不断循环直到在给定的尺度上没有任何对象可以合并为止 ^{[ 15]} . 描述异质性因子间的关系如 图3所示,异质性 h主要由形状因子和光谱因子决定, ω _shape、 ω _color代表形状因子和光谱因子两者间的权重调配,两者之和为1. 形状因子又由光滑度因子、紧致度因子及两者权重总和为1的因子 ω _smooth、 ω _compact计算得到,而光谱因子则由图层权重、图层的标准差和图层数决定.

图  2  多时相影像纹理特征

Figure  2.  Texture feature in multi-temporal images

下载: 全尺寸图片 幻灯片

图  3  异质性因子间的关系

Figure  3.  Relationship of heterogeneity factors

下载: 全尺寸图片 幻灯片

分割过程中需要确定3个因子:分割尺度、波段权重及形状因子 ^{[ 5, 16]} ,形状因子决定了对象形状间的差异;波段权重则描述了不同波段的重要性程度;分割尺度从根本上决定了影像对象的大小. Woodcock等 ^{[ 8]} 认为影像的空间结构可以表现为区域目标大小与空间分辨率之间的关系,并计算影像的局部方差来确定最优分辨率. Drgut等 ^{[ 7]} 将该思想运用到多尺度分割中,通过计算每个层次上对象标准差的均值,引入了方差变化率来描述层次间的变化,并开发了估算尺度参数工具(estimation of scale parameter,ESP),当变化率达到峰值时,认为是该影像的最优分割尺度,计算公式为

ROC =[ $\begin{array}{c}\frac{\mathit{Scal}{e}_l-\mathit{Scal}{e}_{l-1}}{\mathit{Scal}{e}_{l-1}}\end{array}$ ] ×100 (1)

式中:ROC为方差变化率;Scale _l 为目标层次的局部方差;Scale _{l- 1}表示相邻较低层次的方差. 由于图像的复杂性,该方法计算出ROC并不唯一,即存在多个最优尺度 ^{[ 9]} ,需要根据所研究的地物特性进行进一步的确定.

2.3 面向对象的不同地物特征筛选

分割完后的影像由多个多边形对象所组成,依据不同对象的特征差异进行分类. 纹理特征主要使用波段的灰度共生矩阵(gray level co-occurrence matrix,GLCM)及灰度差向量(gray level difference vector,GLDV)的8个常用特征来描述,光谱特征由均值(mean)、最大差分(max diff)、标准差(standard deviation)和亮度值(brightness)等特征来描述,形状特征由长宽比(Length/Width)描述,同时构建了比值植被指数(ratio vegetation index,RVI)、归一化植被指数(normalized difference vegetation index,NDVI)、差值植被指数(difference vegetation index,DVI)和土壤调整植被指数(soil-adjusted vegetation index,SAVI)作为特征加入到分类中 ^{[ 15, 17]} . 纹理特征、光谱特征、形状特征及通过计算得到的各类植被指数,特征众多,如此庞大的特征体系如果全部输入分类器中,必然造成分类效率的下降,而分类结果也并不一定理想,且特征间相关性大、冗余度高,因此需对特征做进一步的筛选. 本文采用特征间的简单相关和多重相关筛选特征,利用相关系数首先剔除相关性较大且方差较小的特征,再通过多重相关性进行筛选. 具体过程如下 ^{[ 18]} .

将各个对象用 n维变量表示,对所有特征进行归一化处理以消除量纲的影响, N个特征的总体协方差矩阵 S可表示为

S= $\begin{array}{c}\left(\begin{array}{cccc}{\sigma }_{11}^2& {\sigma }_{12}^2& \dots & {\sigma }_{1n}^2\\ {\sigma }_{21}^2& {\sigma }_{22}^2& \dots & {\sigma }_{2n}^2\\ ⋮& ⋮& & ⋮\\ {\sigma }_{n1}^2& {\sigma }_{n2}^2& \dots & {\sigma }_{\mathit{nn}}^2\end{array}\right)\end{array}\begin{array}{c}{\sigma }_{\mathit{ij}}^2\end{array}$ = $\begin{array}{c}\frac{1}{n}\end{array}\begin{array}{c}\stackrel{n}{\sum _{k=1}}\end{array}$ ( x _ik-m _i )( x _ik-m _j ) (2)

式中 σ _ij 为第 i个对象特征与第 j个对象特征的协方差. 第 i个和第 j个特征的相关系数为cof _ij = $\begin{array}{c}\frac{{\sigma }_{\mathit{ij}}^2}{\sqrt[]{{\sigma }_{\mathit{ii}}^2{\sigma }_{\mathit{jj}}^2}}\end{array}$ ,当cof _ij ≥0.9时,认为第 i个对象特征和第 j个特征之间存在高度相关性. 通常认为方差大的代表其承载信息量较大,因此在两者之间选取方差较大因子,并利用方差膨胀因子(variance inflation factor, VIF)来检测特征间的多重相关性,将VIF大于10的因子去除,最终得到所筛选的特征,计算方法为

VIF _i= (1 - $\begin{array}{c}{r}_i^2\end{array}$ ) ^{- 1}(3)

式中:VIF _i 为第 i个对象特征的方差膨胀因子; r _i 为以第 i个对象特征为因变量、其他对象特征为自变量的复相关系数.

2.4 遥感影像分类方法

在影像分割和特征选择的基础上,研究选用最邻近分类法 ^{[ 19]} 、支持向量机 ^{[ 20]} 和最大似然分类 ^{[ 21]} 3种分类方法开展人工林地树种分类,对分类结果进行对比分析,并利用混淆矩阵评价分类精度. 考虑到地理实体格局对尺度的依赖性 ^{[ 22]} ,面向对象分类方法依据选取的研究区域特性,在非植被(第1层)和植被(第2层)2个层次上进行信息提取.

3 结果与讨论

3.1 影像分割

研究中首先以较大步长值50进行快速分割,目视观察各地类适合分割的尺度范围. 尺度小于100时,地物已被分割得十分破碎,而尺度大于300时,地类边界没有得到很好地区分,不适用于进行类别的分割提取. 考虑到区域为人工种植林,同种树种聚集程度较高,而不同树种之间界限较为明显,形状相对紧密、规整,紧致度因子和形状因子不应设置得过小. 研究主要针对植被进行提取,由 图1多时相影像地物光谱曲线可以看出,多时相下6月27日影像的近红外上波段植被区分度较高,因此将近红外波段的权重值设置为1.5,其余波段为1.0.

图  4  单时相与多时相不同的分割效果

Figure  4.  Segmetation of single image and multi-temproal images

下载: 全尺寸图片 幻灯片

设定好所需参数后,利用ESP计算最佳分割尺度,当形状因子、紧致度因子均为0.5时,局部方差和方差变化率随尺度增加的变化情况如 图4所示,可看出不同尺度区间上存在多个方差变化率的峰值,对应的最优尺度也不唯一,对比每一尺度参数下的分割效果,最终确定的分割参数如 表1所示.

图  5  最优尺度的计算

Figure  5.  Calculation of optimal segmentation scale

下载: 全尺寸图片 幻灯片

表  1  分割参数表

Table  1.  Parameters of multi-scale segmentation

时相		形状因子	紧致度因子	最优尺度
单时相	第1层	0.4	0.6	220
	第2层	0.3	0.5	164
多时相	第1层	0.4	0.5	170
	第2层	0.3	0.5	114

下载: 导出CSV 

| 显示表格

表1中可以看出,将其他时相影像作为附加波段参与分割过程后,相同参数下多时相影像计算出的最优分割尺度与单时相下的并不相同. 为对比单时相和多时相下的分割效果,设定相同尺度、相同参数因子,对影像进行分割,如 图5所示. 可以看出多时相影像分割的对象更好地描述了不同树种的边界外围,单时相影像异质性较强,分割结果显得更为破碎,生成的对象较多,直接造成后续运算效率的降低. 由于所选择的多时相影像间隔时间仅为几个月,树木经历从3月萌芽到6月底长势旺盛阶段,在此期间认为不同树种边界范围保持不变,即反映在影像上其轮廓位置一定. 而多时相影像分割过程中,根据植被因物候不同引起的光谱、形状信息的差异,可最大程度上保持树种(特别是针叶树种)边界范围的完整性.

3.2 影像特征筛选

文中选取纹理特征、光谱特征、形状特征及通过计算得到的各类植被指数,共构建了60个特征. 在已分割的图像中选择包含各个地类的100个典型对象样本,其中每类地物样本不少于5个,房屋、道路、水体、阴影等较为明显地物样本所占样本总比例为29%,其余植被样本比例为71%. 将在eCognition中构建的所有特征属性值全部导出到SPSS17.0,利用简单相关系数进行初步筛选,探测特征间相关性,剔除相关系数大于0.9且因子方差较小的特征. 再采用VIF对多重相关性进行筛选,将VIF大于10的特征因子去除,单时相、多时相下最终筛选出来的特征如 表2、3所示. 可以看出,特征个数较原来减少许多,对于波段而言,特征波段间自相关性较大,且以第4波段独立性最好,建立的植被指数最终仅保留了方差较大的NDVI. 多时相纹理特征在保留特征中占有较大比例,表明纹理特征与其他特征相关性、复相关性较小,是分类过程中不可或缺的一个重要指标特征.

表  2  单时相影像筛选的特征

Table  2.  Feature selected with single image

最终筛选特征	波段
Length/Width
Mean Layer	NIR
Standard deviation	NIR
GLCM Homogeneity	G
GLCM Contrast	B、NIR
GLCM Entropy	NIR
GLCM Ang. 2nd moment	G、NIR
NDVI

下载: 导出CSV 

| 显示表格

表  3  多时相筛选出的特征

Table  3.  Feature selected with multi-temproal images

最终筛选的特征
Mean_627_NIR	Mean_423_NIR
Mean_317_NIR	GLCM_Mean_423_G
Length/Width	GLCM_Cont_627_NIR
GLCM_Homo_518_NIR	GLDV_Ent_518_R
GLDV_Ent_423_G	GLCM_Cor_423_NIR
NDVI_627	Stddev_317_NIR
GLCM_Ang_423_NIR	GLCM_Mean_627_NIR
GLDV_Mean_317_NIR	Std_dev_627_NIR
注:1)表中数字代表影像日期,R、G、NIR分别代表红、绿、近红外波段;2)表中纹理特征部分采用缩写,分别为同质度(homogeneity,Homo)、对比度(contrast,Cont)、角二阶矩(angular second moment,Ang)、,熵(entropy,Ent)、非相似度(dissimilarity,Dis)、均值(mean)、标准差(stddev)、相关(correlation,Cor).

下载: 导出CSV 

| 显示表格

3.3 树种分类结果

参照外业所采集的样点及目视解译,选择合适、典型的样本数量为210个,每类地物样本不少于15个,使其尽量均匀分布在研究区域内,并按照3:1比例分配训练样本和验证样本,4种分类方法最终得到的结果如 图6和 表4所示(其中金叶榆、普通榆树合并为榆树类别). 从 图6可以明显看出SVM分类、MLC分类中道路与房屋有错分现象,单时相KNN分类、多时相KNN分类图中树木有着较好的聚集范围及清晰的外围轮廓,而对于基于像元分类结果来看,同类像元聚集程度不高,MLC分类效果极为破碎.

仅依靠单时相光谱信息来区分植被的SVM和MLC分类精度分别为49%、43%,Kappa系数分别为0.40、0.33. 而面向对象分类方法的单时相KNN分类总体精度51%,Kappa系数0.42,精度高于传统的监督分类. 参考各树种在不同时相上表现出的光谱、纹理等特征,进行多影像叠加,多时相KNN分类总体精度达到64%,Kappa系数为0.58.

图  6  MLC、SVM、单时相KNN及多时相KNN的分类结果

Figure  6.  Classification of MLC, SVM, KNN in single image and KNN in multi-temproal images

下载: 全尺寸图片 幻灯片

表  4  单时相KNN、多时相KNN、SVM及MLC的分类精度

Table  4.  Accuracy of classification of KNN in single image, KNN in multi-temproal images, SVM and MLC

类别	单时相KNN		多时相KNN		SVM		MLC
	生产者精度	用户精度	生产者精度	用户精度	生产者精度	用户精度	生产者精度	用户精度
房屋	0.64	0.70	0.95	0.90	0.57	0.68	0.53	0.67
道路	0.67	0.65	0.77	0.94	0.66	0.61	0.69	0.59
阴影	0.59	0.95	0.91	0.98	0.50	0.78	0.42	0.81
水体	0.88	0.96	1.00	1.00	0.46	0.82	0.41	0.85
针叶	0.16	0.15	0.65	0.57	0.20	0.16	0.14	0.15
杨树	0.56	0.58	0.69	0.68	0.43	0.44	0.43	0.44
柳树	0.64	0.54	0.67	0.62	0.57	0.51	0.53	0.48
槐树	0.52	0.50	0.60	0.57	0.41	0.38	0.43	0.39
榆树	0.37	0.43	0.73	0.58	0.32	0.40	0.26	0.38
白蜡	0.43	0.54	0.61	0.56	0.42	0.46	0.37	0.44
其他	0.58	0.49	0.46	0.59	0.52	0.42	0.55	0.43
总体精度	51%		64%		49%		43%
kappa系数	0.42		0.58		0.40		0.33

下载: 导出CSV 

| 显示表格

验证结果表明,无论在单时相还是多时相下,房屋、道路、阴影和水体提取精度均优于其他几种地物类型,主要因为此4类地物光谱、纹理特征与其他地物有明显的区别. SVM、MLC分类过程中道路与房屋、阴影与水体有些混淆,而面向对象分类过程中,影像经过分割后形成聚集的一个个对象,在一定程度上保证了房屋、道路、阴影及水体边界的完整性. 另外由于初春、初夏时节,树木还未完全茂盛,阴影较少,盛夏被遮蔽的房屋、道路、水体均有所显露,基于多时相影像分类很好地区分了被掩盖的地物与阴影,面向对象KNN分类对4类地物的提取精度明显提高.

综合来看,单时相基于像元分类方法中,SVM提取的地物精度较好,分类器性能优于MLC. 与基于像元分类相比,面向对象分类方法有着较好的适用性,每类地物信息的提取精度都有所提高,而面向对象方法中多时相影像的加入相对于仅用单时相影像精度也有较好的提升. 研究时间跨度仅仅几个月时间,每类树种分布范围是一定的,即地物的边界轮廓保持不变或变化极少,因此在面向对象分类过程中,多时相影像不仅提供了丰富的纹理、物候信息,而且在影像分割过程提供了较好的边界信息,这对分割后对象边界范围的准确性起到了一定的保证作用,也为后续分类打下了良好的基础,是分类精度提高的另一个重要原因.

树种分类中针叶树种分类精度最低,用户精度、生产者精度均不足0.4,究其原因可能是由于研究区域内针叶树林龄低、树冠小,反映到图像上光谱信息较弱,与区域中灌木较为相似,且与其他树种混合种植更是增加了提取难度. 但多时相物候信息的加入,较大程度上提高了针叶树种的分类精度.

遥感影像上槐树和普通榆树光谱、纹理特征很相似,不易分辨,故两者在所有树种中提取精度较低. 杨树、柳树长势旺盛,上层树冠较为饱满,而白蜡树龄较小,林冠呈现小团状,冠幅差异是造成杨树、柳树区分度较高而白蜡较低的原因.

平原造林中种植林木由于考虑到景观效应,种类繁多、分布复杂,林下地被及分布情况先验知识不足,这影响了树种的分类精度. 虽然在多光谱遥感影像中不同树种幼小林木区分难度较大,但不同树种幼龄林在较窄的光谱波段下存在差异,未来可考虑利用高光谱遥感 ^{[ 23]} 进行幼林树种分类,进一步提高分类精度.

4 结论

本文针对树种混杂、景观差别大、郁闭度低的人工造林地,利用多时相、高空间分辨率卫星影像,通过计算方差变化率,结合目视解译,获取面向对象分类过程中影像的最佳分割尺度;筛选构建的特征,进行了面向对象的多时相影像和单时相影像分类,同时对比了基于像元的分类精度,具体结论如下:

1) 分割过程中通过设定均质性因子,利用ESP计算方差变化率结合目视解译,可更快速地获取适合该影像的最佳分割尺度,且因植被物候信息不同,基于多时相影像树种分割效果优于单时相.

2) 利用简单相关系数进行特征的初步筛选,剔除相关性较大且方差较小的特征,再采用方差膨胀因子对多重相关性进行筛选,最终保留相关性较小同时包含信息量较大的特征用以后续分类,有效减少了特征数量.

3) 基于面向对象下的KNN分类方法,单时相和多时相分类精度分别达到51%、64%,提取精度优于仅利用光谱信息来区分植被的SVM(49%)和MLC(43%)方法,参考各树种在不同时相上表现出的光谱、纹理等特征,多时相影像的加入有效提高了分类精度.