姜萍,赖双双,许婷婷
(1.新疆维吾尔自治区气象服务中心,乌鲁木齐 830002;2. 兰州大学 资源环境学院,兰州 730000;3.广西师范学院 国土资源与测绘学院,南宁 530001;4.新疆维吾尔自治区气象台,乌鲁木齐 830002)
近年来,我国针对积雪面积信息提取的遥感研究主要围绕波段数较多、光谱信息相对比较全面的Landsat TM/ETM+,SPOT,MODIS等卫星数据开展[1],并且具备全天时、全天候以及穿透探测能力的SAR卫星影像也逐渐参与其中,在积雪深度、雪水当量反演方面优势突出而获得广泛应用与研究[2-4]。与此同时,针对上述数据的多种雪盖信息提取方法得到迅速发展。张永宏等[5]利用归一化雪盖指数法(normalized difference snow index,NDSI)和多波段综合阈值实现积雪判别,研究出了适用于FY-3/VIRR资料的积雪判识方法,通过实例分析验证算法有效可行。侯慧姝等[6]利用NDSI算法和综合阈值判别法对内蒙古自治区大范围降雪进行积雪信息提取。Gu等[7]基于图像匹配技术,采用NDSI、NDVI及多波段综合阈值法,提出了一种综合MODIS和AMSR-E数据的积雪提取方法。瞿娟等[8]利用MODIS数据,提出了结合混合光谱分解的积雪分量及灰度共生矩阵提取的纹理特征的支持向量机(support vector machines,SVM)分类方法,对研究区积雪面积信息提取进行了研究。Tang等[9]通过修正地表反射率、提出归一化云指数(normalized difference cloud index,NDCI)区分云/雪像元、引入MOD11_L2地表温度产品进一步提高方法精度等方式,提出了一种改进的基于MODIS影像的喜马拉雅地区雪盖提取方法。何咏琪等[10]和蒋友严等[11]分别以Landsat ETM+与MODIS提取的积雪面积为标准,分析确定了适用于HJ-1B卫星数据的NDSI阈值,并结合SNOMAP算法提取积雪面积,验证了算法的可靠性。郭金金等[12]对全极化SAR数据进行目标分解提取积雪极化特征,再利用J-M距离进行特征选择,分析不同极化特征对积雪的可分性,最后利用最优特征集和SVM进行积雪识别,能够获得较好实验效果。众多方法中,阈值法及基于反射率特性的NDSI法在积雪盖度信息遥感提取中应用最为广泛。
2012年1月9日,我国资源三号(简称ZY-3)卫星成功发射,较之上述卫星数据,其获取的4波段(其中1~3波段为可见光波段(R、G、B)、第4波段为NIR波段)多光谱影像具有更高的空间分辨率[13]。如何针对国产ZY-3卫星多光谱影像空间分辨率高但光谱范围较窄、缺少NDSI计算中所需短波红外波段这一光谱特性,寻求一种快速、准确的雪盖信息提取方法极具研究意义。本文基于积雪与非积雪地类在ZY-3卫星多光谱影像中光谱反射特性的差异性,以及各地物的NDVI特性,提出一种综合NIR波段与NDVI的决策树雪盖信息提取方法。实验表明,该方法能保证雪盖信息提取的完整性和准确性,在雪盖信息提取方面具有一定的优越性。
本文实验采用2012年2月2日获取的ZY-3 卫星高分辨率多光谱影像作为数据源,分辨率5.8 m,云量6%,成像区域为新疆乌鲁木齐市北部与昌吉回族自治州接壤处,范围为87°37′57″E~87°47′45″E,44°7′57″N~44°19′27″N。研究区地处北疆地理中心,属于温带大陆性气候,夏季短促而炎热,冬季漫长而严寒,每年都要经历3~5个月的降雪期,积雪资源非常丰富,占全国积雪资源的1/3,便于开展积雪研究且应用价值较大。研究区内土地覆盖类型主要包括积雪、建筑物、道路、裸土和少量薄冰等(基本无云层干扰)。研究区地理位置见图1。右图中黑色矩形框所示区域为本文所用实验区。
图1 ZY-3多光谱影像及其地理位置
ZY-3卫星是我国首颗高精度民用立体测绘卫星,设计寿命为5年,主要任务是长期、连续、稳定和快速地获取全国高分辨率立体影像和多光谱影像。该卫星同时搭载2.1 m分辨率的全色传感器和5.8 m分辨率的多光谱传感器,其影像参数如表1所示。
表1 ZY-3多光谱影像波段参数与绝对定标增益系数
ZY-3卫星多光谱影像的预处理包括波段表观辐亮度和波段表观反射率的计算2个方面。
1)辐射校正。采用绝对辐射定标系数(中国资源卫星应用中心网站获得),将ZY-3卫星影像各传感器通道观测值DN转换为传感器入瞳处表观辐亮度值。
2)Flaash大气校正。根据ZY-3卫星多光谱数据的光谱响应函数,采用Flaash大气校正方法消除大气和光照等因素对地物反射的影响,获得地表反射率。
研究表明,积雪在可见光波段(0.4~0.75 μm)反射率较高,对太阳光谱几乎达到全反射,在0.5 μm附近出现反射峰,反射率高达80%以上[14],如图2所示。随着波长的增加,积雪反射率逐渐降低,但仍在可见光范围内保持较高的反射率,从而在影像中表现为浅白色。不被积雪覆盖的裸土、道路、建筑物等地类,因反射率较低在影像中表现为暗色。
图2 积雪反射率光谱特性曲线
为进一步分析ZY-3卫星多光谱影像上积雪与非积雪地物类型的光谱差异性,分别在实验区1中选择积雪样本,裸土、道路、建筑物、薄冰等非积雪样本各400个,对其在4个波段中的地表反射率进行统计,结果如图3所示。
图3 ZY-3卫星影像中不同地类表观反射率曲线
图3中显示,积雪覆盖程度较大的区域(深雪),表现为强反射特性,其表观反射率远大于非积雪地类。而积雪覆盖程度较小的区域(薄雪),表观反射率有所下降,并在可见光波段与薄冰相近。其他4类非积雪地类在4个波段中的表观反射率均远小于积雪覆盖地类。其中,第四波段(即NIR波段)中,无论是深雪区,还是薄雪区,其表观反射率都明显高于其他地类。由此可见,ZY-3卫星多光谱影像NIR波段具有识别积雪信息的能力。
NDVI是监测地区或全球植被和生态环境变化的有效指标,由NIR波段和红波段反射率比值经非线性归一化处理得到。
NDVI取值在[-1,1]范围内。在植被覆盖下,NDVI>0,且随植被覆盖度增大而增大,因此该指数被广泛应用于植被提取中。而研究表明,对于除植被以外的其他陆地地表覆盖类型,由于其在可见光波段、NIR波段上反射率差异的不同,使得NDVI同样具有一定的指示作用。如岩石、裸土等在两波段上反射率相近,NDVI≈0;积雪、云、水、阴影和薄冰等在可见光波段反射率高于NIR波段,NDVI<0[15]。有积雪覆盖的地表,NDVI<0即为积雪的NDVI特性,目前已有学者将此特性用于雪盖制图方法研究,并证实了方法的可行性[16]。
本文针对国产ZY-3影像分辨率高但光谱范围较窄,缺少NDSI计算中所需短波红外波段这一特点,以对ZY-3卫星多光谱影像中积雪与非积雪光谱反射特征分析、积雪的NDVI特性分析为基础,提出了一种综合NIR波段与NDVI的决策树雪盖信息提取方法。
方法首先充分利用积雪在NIR波段上反射率远高于其他地类的特性,选取合适的阈值T,对ZY-3影像进行二值化处理。阈值T选取过程中,考虑到大量建筑物楼顶亦被积雪覆盖,为使该部分雪盖信息被最大程度地保留,故应选取稍小的阈值。参考图3中统计结果,选择T=0.74作为NIR二值化的阈值,进而获得雪盖信息初步提取结果;其次,为进一步去除初步提取结果中含有的薄冰、未覆雪建筑物等非积雪信息,以NDVI<0为判别指标对该结果进行二次阈值化,即充分利用积雪的NDVI特性,实现最佳雪盖信息的提取。图4为本文决策树雪盖信息提取的技术流程。
图4 决策树雪盖信息提取技术流程
图5为实验区1对应的431波段假彩色合成影像、最大似然监督分类方法和本文方法提取的雪盖信息对比。其中,人为参与样本选择的最大似然监督分类结果将作为参考数据,对本文方法的提取效果进行定性与定量分析。
图5 实验区1雪盖信息提取对比
从图5(a)、图5(c)对比中可直观看出,整体来看本文提出的提取方法能够较好地将实验区内积雪像元提取出来。实验区1中积雪覆盖程度较大的青格达湖、大面积农田、裸地、及部分屋顶覆雪的建筑物等在提取结果中获得了较好的保留;而实验区1左上方的五家渠市、右侧的羊毛工等村镇、左下方青格达湖景区一侧的山林等大量未覆雪地类也得到了较为完整地表现。由图5(b)、图5(c)对比可见,本文方法提取效果与最大似然分类法吻合程度较高。
采用总体精度和Kappa相关系数来进一步定量评价本文方法提取效果。表2为本文方法与最大似然分类法提取结果的混淆矩阵,可见,本文方法提取结果总体精度为94.97%,Kappa系数为0.884,表明方法能较准确而完整地识别出实验区内的积雪像元。综合图5及表2可见,本文方法能保证雪盖信息提取结果的完整性和准确性,在雪盖信息提取方面具有一定的优越性。
表2 本文方法提取结果与最大似然分类法提取结果混淆矩阵
选取实验区1中建筑物密集区域、含阴影和裸土区域、薄冰以及积雪覆盖程度较小的农田4个区域进行提取结果细节对比,具体比较如表3、表4所示。
从表3、表4中可以看出,在建筑物、阴影、裸土和薄冰等区域本文方法与作为参考数据的最大似然分类法提取结果一致性高。对于区域1、区域2,本文方法能够有效识别和划分屋顶有无积雪覆盖的建筑物、城市内不同等级的道路、高大楼层的阴影以及冬季草木枯萎而裸露的城市绿地等地物类型,同时清晰地提取出城市中大量人工结构的轮廓等细节信息,漏提及误提现象较少,提取总体精度分别为93.56%、93.68%,Kappa系数均在0.85左右。区域3中,本文方法较最大似然分类法提取的雪盖面积有所增加,且增加部分集中在薄冰区边缘一带。经分析,其原因在于薄冰区在NIR波段的光谱反射率基本处于判别阈值0.74之上,而冰边缘通常为积雪所覆盖,致使其对应NDVI<0,从而被划分为积雪一类。结合现实考虑,认为该部分误分可视为合理性的划分。同时,因该部分像元个数所占比例较小,对提取结果总体精度影响不大。区域4中对于被薄雪覆盖的农田存在较严重的漏提,总体精度仅有87.71%,Kappa系数0.687。经分析,NDVI>0是引起该现象的主要原因。如何有效解决该部分漏提问题,有待进一步研究。
为进一步验证本文所选分割阈值T的代表性与普适性,选取昌吉市边缘上玉堂村周边地区作为实验区2,直接采用本文算法进行覆雪面积的提取。由结果(图6(b))可见,实验区2中主要覆雪地类提取效果与真实情况吻合程度较高,而主要道路、城乡建筑、和部分裸土均被正确判别为非覆雪区。由此证明,在保持分割阈值不变的情况下,本文算法依然能够取得较好效果。
表3 雪盖提取结果局部对比
表4 局部区域提取结果混淆矩阵
图6 实验区2雪盖信息提取对比
本文以ZY-3卫星多光谱影像为数据源,综合NIR波段与NDVI提出了决策树雪盖信息提取方法。通过对提取结果进行评价分析发现,本文方法能较准确有效地提取出研究区中的雪盖信息,总体精度可达94.97%。由局部区域对比发现,本文方法受建筑物、阴影和裸土的影响较小,能在保留各地类细节信息的同时获得较高的雪盖信息提取精度;在薄冰区边缘,由于积雪覆盖程度增大,其对应NDVI<0引起提取结果中该部分雪盖面积合理性的增加,但总体影响较小;被薄雪覆盖的农田在本文方法中存在较为严重的漏提现象。今后的研究重点将集中于消除现有漏分、误分问题,提高阈值选取的精确性,同时开展对含云量较多的影像的积雪识别研究,为ZY-3卫星多光谱影像在积雪资源方面的研究和应用提供支持和参考。