字 李, 谢福明, 舒清态, 吴 荣
(西南林业大学林学院,云南 昆明 650224)
与传统的遥感数据相比,高光谱数据具有分辨率高、波段多、信息含量丰富的独特优势,然而其庞大的数据量也给研究和应用带来了困难.为此,国内外学者做了大量关于高光谱树种精细判别的研究.王志辉[1]、洪娇[2]等对比了不同树种的原始光谱、一阶微分光谱、二阶微分光谱曲线图,发现区分不同树种的波段主要位于近红外波段,并且利用欧式距离法对波段选择结果进行检验. 褚西鹏等[3]采用平滑处理提高了叶片高光谱数据树种分类的精度.丁丽霞等[4]利用包络线去除法对高光谱数据降维处理,有效解决了高光谱数据冗余的问题.刘秀英等[5]发现原始数据对数变换的一阶微分的杉木和马尾松的判别精度最高.相比于国内,国外对于树种精细判别的研究相对较早,应用更为广泛.Prospere et al[6]在判别牙买加热带湿地树种的研究中发现弹性网规则化广义线性模型更适合于基于特征选择的光谱指数. Meng和Dennison[7]使用特征选择技术,采用决策树算法分类比较冠层尺度的高光谱数据,减少数据冗余和计算时间,提高了树种分类精度.但目前利用非成像高光谱遥感数据精细判别森林针叶树种的研究较少,对生态脆弱区典型针叶树种的判别就更少.
本试验以香格里拉市生态脆弱区云南松(Pinusyunnanensis)、高山松(Pinusdensata)、云杉(Piceaasperata)、冷杉(Abiesfabri)4种典型针叶树种为研究对象,在利用ASD Field Spec 3地物光谱仪野外采集光谱数据和预处理的基础上,对实测的叶片高光谱数据进行原始光谱、一阶微分光谱、二阶微分光谱分析,借助Fisher判别分析方法对4种典型针叶树最佳波段窗口和判别精度进行验证,旨在为中大尺度机载、星载高光谱遥感树种精细分类提供指导.
位于滇西北的香格里拉市地处世界自然遗产“三江并流”核心区(99°20′—100°19′E、26°52′—28°52′N),平均海拔3 459 m.主要气候类型为山地寒温带季风气候,降雨量充足,优越的水热条件为当地森林的生长提供了良好的环境.香格里拉市林地面积9 509 km2,森林覆盖率74.99%,主要有高山栎、云冷杉、落叶松、云南松、高山松、华山松等10种植被类型.本试验所选的4种树种均为当地典型针叶树种,合计占林地总面积的60.20%,阴坡、阳坡均有分布.其中云南松林面积984 km2,占林地总面积的10.35%,分布海拔为1 700~2 800 m,主要土壤为褐土、红壤、黄棕壤、亚高山和高山草甸土,上限与亚高山的高山松林犬牙交错;高山松林面积1 742 km2,占林地总面积的18.32%,分布于云杉、冷杉林下限,分布海拔为2 800~3 500 m,主要土壤为棕壤、亚高山和高山草甸土;云杉林面积956 km2,占林地总面积的10.05%,分布于冷杉林下限、高山松上限,分布海拔为2 900~3 500 m,主要土壤为亚高山和高山草甸土;冷杉林面积2 044 km2,占林地总面积的21.50%,分布上限接高山杜鹃灌木林或高山草甸,下限接云杉林,分布海拔为3 000~4 500 m,主要土壤为暗棕壤、棕色针叶林土、亚高山和高山草甸土.
本试验采用ASD Field Spec 3地物光谱仪采集4种针叶树种的叶片光谱反射率,该设备在350~1 500 nm波长范围内分辨率为3 nm,采样间隔为1.4 nm;在1 500~2 500 nm波长范围内设备分辨率为10 nm,采样间隔为2 nm.
由于测定的树种比较高大,故采用离体测量的方法采集数据[8],即先用高枝剪剪下向阳面健康的树枝,将松针紧密地铺在黑布上,再进行数据采集.采集时间为2015年11月18—21日11:00—16:00,该时段天气晴朗少云无风.选用10°光谱仪视场角,测定时探头垂直向下,与采集样本的垂直距离以松针全部落入视场范围内为准.每次采集光谱数据前,都进行一次标准的白板校正,并根据当地的天气变化实时优化系统,取10个光谱值的平均值作为一个采样光谱样本.每个树种选取10棵树,从每棵树的3个高度(2、4、6 m)采集6个向阳面的健康树枝.以枝条为单位,每种树采集60个样本,共计240个样本数据.
1.3.1 异常光谱数据处理 在处理异常光谱时,首先采用View SpecPro 6.0光谱处理软件导出由美国ASD Field Spec 3地物光谱仪测量的原始光谱反射率曲线,每次显示10条光谱反射率曲线,删除与整体样本偏离较大的样本和异常光谱样本,重复以上步骤直到整个样本数据分布均匀.然后检查样本数据的所有波段(350~2 500 nm),删除受空气中的水汽以及仪器的系统误差影响较大的波段(1 350~1 400、1 800~1 950、2 300~2 500 nm)的数据,保留质量较高的波段(350~1 349、1 401~1 799、1 951~2 299 nm)的数据.
1.3.2 均值处理 由于不同针叶林类型不同高度叶片的生化参数空间分布不均匀,导致各高度对应的光谱反射率也有所变化.将不同高度枝条测得的光谱反射率值取平均作为一棵树的光谱反射率值,某一树种的光谱反射率值是所有树种样本的光谱反射率平均值.
1.3.3 平滑去噪处理 研究表明,采用Savitzky-Golay滤波算法对数据进行平滑去噪,不仅能有效去除脉冲噪声和白噪声,还可以保留原始高光谱数据的奇异特征并对移动窗口的曲线进行拟合,从而提高光谱数据特征波段提取的精度[9],计算公式如下:
(1)
式中,Yi为滤波后的光谱曲线,Yi+1为原始光谱曲线,Cj为滤波系数,N为滑动窗口所包括的数据点2m+1.对比分析本试验光谱曲线实际情况,滑动窗口N取值为25.
1.3.4 数据统计分析 利用Excel 2013软件对预处理得到的原始光谱反射率、反射率一阶微分和反射率二阶微分结果逐个进行分类整理,得到3张相对应的样本光谱特征表.随机选取2/3(云南松、高山松、云杉和冷杉各40个)样本进行敏感波段分析,选取1/3样本(云南松、高山松、云杉和冷杉各20个)进行SPSS判别分析及精度检验.
由于受仪器自身和大气的影响,野外采集的非成像高光谱数据之间普遍存在系统误差以及大气吸收、散射和辐射等背景噪声.利用光谱微分法不但可有效消除这些误差和背景噪声,还可分辨原始光谱的重叠光谱,增强微分光谱曲线上的细微变化,提取出可判别树种的光谱吸收峰参数,以达到树种的判别精度[10-13].
1.4.1 光谱一阶微分 光谱一阶微分代表原始光谱曲线在各个波段的斜率,反映了原始光谱反射率升降的程度,光谱反射率升降分别用正值和负值表示[11],计算公式如下:
(2)
式中,λt为第波段,FDRλt为波段t和波段t+1之间的光谱一阶微分,Rλt、Rλt+1为相对应的波段t、t+1处原始光谱反射率值,Δλ为1.
1.4.2 光谱二阶微分 光谱二阶微分代表了原始光谱曲线的弯曲状况,原始光谱曲线凹凸分别用用正值和负值表示,相比于光谱一阶微分,光谱二阶微分更加有效地扩大了光谱间细微的差别[14],计算公式如下:
(3)
式中,λt为第波段,SDRλt为波段t和波段t+2之间的光谱二阶微分,Rλt、Rλt+1、Rλt+2为相对应的波段t、t+1处原始光谱反射率值,Δλ为1.
Fisher判别分析的基本思想是投影,将k组p维数据投影到某一个方向,使得它们的投影对组与组之间来说差异尽可能地大[15].给定一个样本,要判断它属于哪个总体,其判别函数为:
(4)
实际应用中,要从贡献率最大的函数选起,依次选择第一个函数,第二个函数,第三个函数……使函数的累积贡献率达到一定的水平.利用Fisher判别函数判别后,使各组间的差异尽量大,而组内各个样本的差异尽量小,组间差异和组内差异的比值越大越好.Fisher判别分析可以筛选出树种判别能力最强的波段并显示判别精度,对多波段高光谱数据的判别效果更佳[16].
本试验对香格里拉4种典型针叶树种进行判别时,先基于原始光谱曲线、光谱一阶微分曲线和光谱二阶微分曲线选出特征波段,分别对应到原始光谱的80个样本,筛选出每组特征波段所对应的原始光谱反射率值,利用SPSS 17.0分别对3个特征组处理结果进行Fisher判别分析.Fisher判别分析时分组变量选择“树种”,自变量为对应的3种数据转换方法所选出来的特征波段;采用步进式方法,使用马氏距离进行逐步变量选择,选用F值作为逐步判别停止的标准, 进入变量和剔除变量的F值分别为3.84和2.71.由于4种树种的所选样本数量相等,在进行逐步判别时,每种树种的权重值相等,在进行4种典型针叶树种判别的基础上,针对根据空间分布特征分开的两组树种进一步判别.
图1 4种典型针叶树种原始光谱曲线Fig.1 The original spectrum curve of 4 tree species
2.1.1 4种典型针叶树的光谱曲线 4种典型针叶树种的原始光谱曲线如图1,整体上4种树种的曲线状态较为相似,符合绿色植被的光谱曲线.由于叶绿素强烈的吸收作用和反射作用,在490、550、680 nm附近呈现“蓝谷现象”、“绿峰现象”、“红谷现象”,在710~800 nm波段反射率陡峭上升,在800 nm附近呈现“红边现象”.但在特定的波段仍然存在差异,如在“绿峰”附近,冷杉的反射率明显低于其它树种;在“红边”附近,云南松的反射率明显高于其它树种.
由光谱一阶微分曲线(图2)可知,在730 nm附近原始光谱反射率变化程度最大,且云杉的变化程度明显高于其它三种树种;在1 140 nm附近,其反射率变化程度也存在明显差异.由光谱二阶微分曲线(图3)可得,原始光谱曲线凹凸比较明显的波段主要位于690、740、1 160、1 410 nm附近.
图2 4种典型针叶树种光谱一阶微分线Fig.2 The first derivative spectrum of 4 tree species
用于树种判别,以10 nm为步长,在原始光谱曲线、一阶微分光谱曲线、二阶微分光谱曲线分别选出21、11、9个波段,具体波段详见表1.
表1 4种典型针叶树种不同光谱处理方法的显著差异波段Table 1 Bands which showed significant difference by 3 spectra processing methods
2.1.2 云南松、高山松的光谱曲线 在海拔2 800 m附近,云南松、高山松交替分布,二者具有极其相似的叶片形态和空间分布特征,所以对其进行做进一步的判别.原始光谱差异最明显的波段为800~1 120 nm(近红外波段),在730 nm附近二者原始光谱反射率变化程度最大(云南松大于高山松),原始光谱曲线凹凸比较明显的波段位于690 nm和740 nm附近,以10 nm为步长从中选取光谱差异性较大波段用于云南松、高山松的判别(表2).
表2 云南松、高山松不同光谱处理方法显著差异波段Table 2 Bands that showed significant difference by 3 spectra processing methods when distinguishing P.densata from P.yunnanensis
2.1.3 云杉、冷杉的光谱曲线 对比云冷杉的原始光谱曲线(图1)、光谱一阶微分曲线(图2)、光谱二阶微分曲线(图3)可以看出,二者原始光谱在550(绿峰)、1 190、1 650 nm附近差异明显,在730 nm附近二者原始光谱反射率变化程度最大(冷杉大于云杉),原始光谱曲线凹凸比较明显的波段位于690 nm和740 nm附近,以10 nm为步长从中选取差异性较大的波段用于云杉、冷杉的判别(表3).
表3 云杉、冷杉不同光谱处理方法显著差异波段Table 3 Bands that showed significant difference by 3 spectra processing methods when distinguishing P.asperata from A.fabri
2.2.1 4种典型针叶树种的判别 4种针叶树不同光谱处理方法的树种判别精度及入选波段见表4.原始光谱共入选7个波段,其中5个位于近红外波段,1个位于绿光波段,1个位于蓝光波段,最佳波段窗口位于近红外波段(980~989 nm).一阶微分光谱共入选8个判别波段,3个位于近红外波段,1个位于红光波段,2个位于绿光波段,2个位于蓝光波段,最佳波段窗口位于蓝光波段(415~424 nm).二阶微分光谱入选8个波段, 其中5个位于近红外波段,2个位于绿光波段,1个位于蓝光波段,最佳波段窗口位于近红外波段(960~969 nm).3种数据处理方法树种判别的总体精度均不低86.3%,说明利用Fisher判别分析对特征波段的降维效果显著,可用于叶生态脆弱区典型针叶林树种的判别.此外,一阶微分光谱和二阶微分光谱的总体判别精度(分别为97.5%和98.8%)都高于原始光谱的判别精度(86.3%),说明对原始光谱进行微分变换,尤其二阶微分变换可以有效地提高针叶林树种精细判别的精度.
表4 4种典型针叶树种不同光谱处理方法的树种判别精度及入选波段Table 4 Accuracy of different spectrum processing methods and candidate bands for spectral discrimination
2.2.2 云南松、高山松的判别 云南松、高山松判别精度及判别函数所入选波段见表5,原始光谱入选7个波段用于判别,其中6个位于近红外波段,1个位于蓝光波段,最佳波段窗口位于近红外波段(870~879 nm);一阶微分光谱入选的8个波段中,4个位于近红外波段,1个位于红光波段,1个位于绿光波段,2个位于蓝光波段,最佳波段窗口位于近红外波段(1 020~1 029 nm);对于二阶微分光谱入选的5个波段, 1个位于近红外波段,1个位于红光波段,2个位于绿光波段,1个位于蓝光波段,最佳波段窗口位于绿光波段(530~539 nm).相比于4种树种的总体判别精度,其3种数据处理方式总体判别精度明显提高,所选波段可用于两种树种的判别.
表5 云南松、高山松不同光谱处理方法树种判别精度及入选波段Table 5 Accuracy of different spectrum processing methods and candidate bands for spectral discrimination when comparing P.densata and P.yunnanensis
2.2.3 云杉、冷杉的判别 原始光谱共入选5个波段用于判别云杉和冷杉(表6),其中3个位于近红外波段,2个位于绿光波段,最佳波段窗口位于绿光波段(540~549 nm);对于一阶微分光谱,入选的5个波段中,1个位于近红外波段,2个位于红光波段,1个位于绿光波段,1个位于蓝光波段,最佳波段窗口位于绿光波段(520~529 nm);对于二阶微分光谱,入选的7个波段中, 3个位于近红外波段,2个位于红光波段,2个位于绿光波段, 最佳波段窗口位于近红外波段(1 150~1 159 nm).
表6 云冷杉不同光谱处理方法树种识别精度及入选波段Table 6 Accuracy of different spectrum processing methods and candidate bands for spectral discrimination when comparing P.asperata and A.fabri
香格里拉市地处青藏高原南延部分和横断山脉纵谷地带,特殊地貌和寒冷气候导致了脆弱的生态环境,境内森林覆盖率高,树种比较单一.区分不同针叶树种所选择的波段宽度都为10 nm,说明要区分不同的针叶树种,应充分利用高光谱数据波段窄的优势,挖掘出不同针叶树种高光谱信息特征的差异性.已有许多高光谱遥感影像能提供足够窄且覆盖近红外波段的数据,如 MODIS 数据(9 个近红外波段)、 Hyperion数据(216 个近红外波段).
本试验分析了4种不同树种的实测光谱曲线,虽然整体上不同的树种的光谱特性十分相似,符合绿色植被的光谱曲线,但也存在着微小的差异.为有效地对4种针叶树种精细判别,将其分为云南松、高山松和云冷杉两组进一步判别.4种树种、云冷杉和云南松、高山松差异最显著的原始波段分别为980~989、540~549、870~879 nm;利用光谱微分变换处理能够增强这些差异,一阶、二阶微分变换后其差异最显著的波段分别为415~424、520~529、1 020~1 029 nm和960~969、1150~1159、530~539nm;利用光谱微分变换处理能够增强这些差异,提高光谱判断的精确度. 不同谱微分换方法选择的波段有些不同,但用于识别不同树种的波段大都位于近红外波段.Fisher判别方程入选波段中,510~559和680~739 nm是本试验最重要的光谱区域,与任何其它分区相比,510~559 nm包含最多的多次选择波长,这是与叶绿素有关的反射作用影响差异所致;680~739 nm区域包含了多次选择的最高频率的连续波长,这归因于近红外波段的重要性,这与Jones et al[17]、Cross et al[18]等国外研究成果较为一致.虽然实验结果较为理想,但也存在以下不足之处:
(1)在样本数据的采集时,考虑到针叶林类型的树木冠层过高,只是借鉴了前人的研究成果,利用黑布和成簇的松针来模拟树木冠层的环境,而并未采集真实冠层的非成像高光谱数据.在今后的研究之中,若能获取实测的冠层非成像高光谱数据,将更好的为中大尺度机载、星载高光谱遥感树种精细分类提供指导.此外采集的数据由于受到水汽波段的影响剔除了(1 350~1 400、1 800~1 950、2 300~2 500 nm)的数据,而借鉴丁丽霞等[4]的研究结果,选择出的最佳波段中有1 920 nm左右的波段,室外采集的数据一定程度上降低了树种判别的精度.
(2)本试验的数据采集时间为11月份,未能检测到叶片光谱树种判别的最佳时相,基于不同时相的4种树种的叶片光谱的动态变化分析需要进一步探索.而且不同树种具有不同的最佳判别年龄,今后的研究中,可采用树种不同年份的光谱数据,综合各年份的光谱特征进行树种的精细判别.