利用多时相遥感影像分析密云县植被动态变化1)

赵泰安 彭道黎

(北京林业大学，北京，100083)

植被是陆地表面主要的覆盖类型，作为反映区域环境的标志之一，植被是环境重要组成及影响因子，也是人们研究的主要对象[1]。植被盖度指植物群落总体或个体地上部分的垂直投影面积与样方面积之比的百分数。植被盖度测量的传统方法是地面测量，最简单的测量方法就是目估法[2－3]，其缺点是主观性太强，更加客观的测量方法有样方法[4]、样带法[5]、样点法[2，6]等，以及借助采样仪器的测量方法，如空间定量计、移动光量计等[2]。这些方法虽然提高了测量的精度，但野外操作不便，并且成本较高，难以在大范围内快速提取植被盖度。相对于前几种方法，利用遥感手段对地表植被的提取相对简便，一方面在于遥感数据的获取相对容易，成本低廉;另一方面也有利于在大范围内适时的监测植被动态变化。一般来说，利用多时相遥感数据进行变化分析的方法最常用的有两种:逐个像元对比和分类后对比。逐个像元对比是将多时相的遥感图像迭合在一起，对每个像元进行逐一比较，无需进行图像分类，如果像元之间亮度值差别较大，表明该像元发生了的变化。分类后对比，首先对各个图像进行分类，然后将分类结果迭合起来进行对比，从而获得各个像元的变化情况;逐个像元对比由于无需进行分类，因此可以尽量减少分类带来的误差，但从该方法所获得的结果只能知道某一个像元是否发生了变化，无法知道发生变化的内容。分类后对比可以清楚地获得每一个像元具体的变化性质，但是该方法由于要首先进行分类，其结果受分类误差的影响，从而导致普遍出现夸大变化的现象[7]。

由于上述两种变化分析方法各有特色和侧重点，在方法的直观性、可操作性以及定量分析等方面仍有值得改进之处。本文就目前多时相遥感数据进行变化分析的两种常用方法进行一个简要的探讨，即用归一化植被指数(NDVI)逐个像元对比和像元二分法植被覆盖度分类后进行对比，分析密云县植被覆盖度的动态变化，为密云县气候变化和水土保持方面提供科学依据。

1 数据来源及研究区域概况

本文研究选用两幅TM、一幅ETM影像，条带号为123/32，时间分别是1992年9月7日、2009年9月22日和2005年9月3日。数据来源于中国科学院计算机网络信息中心国际科学数据服务平台(http://datamirror.csdb.cn);影像获取的月份一致，有利于进行遥感分析。此外还有密云1∶5万地形图，密云县一类清查资料。

密云县位于北京市东北部，属燕山山地与华北平原交接地。密云县西起东经116°39'33″～117°30'25″，东西长 69 km;南起北纬 40°13'7″～ 40°47'57″，南北宽约64 km。东南至西北依次与北京市的平谷区、顺义区、怀柔区接壤，北部和东部分别与河北省的滦平县、承德、兴隆县毗邻。全县总面积2229.45 km2，占全市面积的13%，是北京市土地面积最大的区县。全县山区面积占79.5%，1771.75 km2;平原面积占11.8%，263.4 km2;水面面积占 8.7%，194.3 km2。属于半湿润大陆性季风气候，冬冷夏暖，年平均气温为10.8℃，年平均降水量为657 mm。

2 研究方法

图像校正最主要是辐射校正、系统粗校正和几何精校正，本次试验所用的数据均经过辐射校正和系统粗校正，只需要对下载的数据进行几何精校正和对ETM图像进行去条带处理。由于landsat－7 ETM+机载扫描行校正器(SLC)故障导致2003年5月31日以后获取的图像出现了数据条带丢失，为了使ETM数据能够使用，必须对其进行去条带处理。在处理中所运用的模型是多影像局部自适应回归分析模型，既利用多景不同时相的遥感数据、采用局部回归分析方法对一景影像进行缝隙填充，回归区域面积为变化值，提取局部区域相关性最大的区域进行回归分析，选择局部区域面积最小，相关性最大的区域进行回归分析。处理结果表明，该模型最最大程度提高了影像修复质量，能满足数据质量的要求。

本研究运用ERDAS IMAGINE9.2中的数据预处理(date prep)模块下的几何校正模块(geometric correction)完成几何校正，借助北京市1∶5万地形图及部分GPS野外实测数据，依据“多点均匀分布、选择相对固定及明显的地物点”。在整幅1992年的TM影像上选取21个控制点，对于不符合精度要求的控制点进行删除、重选，直至符合要求。经RMS检验，误差小于一个像元，采用立方卷积法进行重采样，产生地面分辨率30 m，投影系统为高斯——克吕格投影，北京——54坐标系的几何精校正影像。其余两幅TM影像均以校正好的1992年的TM影像为基准进行影像间的配准，利用影像——影像的方式，达到几何配准的目的。

获取密云县的行政边界图，对其进行投影变换和坐标配准，使其与遥感影像坐标相吻合。首先利用ERDAS的Vector将矢量多边形转换成栅格图像文件，再通过掩膜运算(Mask)实现图像不规则裁剪。本次试验所采用的是ERDAS矢量栅格转换功能和掩膜运算得到密云县的遥感影像。

3 结果与分析

3.1 归一化植被指数的提取

归一化植被指数计算公式如下:

其中，RNI表示近红外波段，R表示红色波段。本文所采用的是Landsat TM数据，因此公式(1)可以表达为:

利用ERDAS遥感处理软件，可以分别从3幅遥感影像中提取出归一化植被指数影像用ERDAS提取3幅归一化植被指数影像的值分别为:－0.7037＜NDVI(1992)＜0.6363，－0.67089＜NDVI(2005)＜0.597，－0.4359＜NDVI(2009)＜0.72093。为了便于对3个时期归一化植被指数合成图进行分析，用ERDAS中Model Maker功能将归一化植被指数图像中像元灰度值归一到0～255之间，得到3幅影像归一化的NDVI值分别为5＜NDVI(1992)＜253，4＜NDVI(2005)＜244，0＜NDVI(2009)＜255(图1)。

3.2 假彩色合成原理

RGB色彩就是常说的三原色，R代表Red(红色)，G代表 Green(绿色)，B代表 Blue(蓝色)。RGB色彩模式使用RGB模型为图像中每一个像素的RGB分量分配一个0～255范围内的强度值。我们将1992年、2005年及2009年3幅单波段影像进行RGB假彩色合成，利用颜色的变化得出密云县植被的动态变化情况(图2)。

3.3 植被盖度遥感估算像元二分模型[8]

像元二分模型最大的优点就是削弱了大气、土壤背景与植被类型等的影响。在前人研究的基础上可知，像元二分模型适用于密云县植被盖度的研究。假设一个像元的信息可以分为土壤与植被两部分，通过遥感传感器所观测到的信息S，就可以由植被信息和土壤信息两部分组成，这里我们将植被信息表示为Sv，土壤信息为SS，那么像元的信息S可以分解为两部分。

对于一个由土壤与植被两部分组成的混合像元，像元中由植被覆盖的面积比例，即为该像元的植被盖度FC，而土壤覆盖的面积比例为1－FC。设全由植被所覆盖的纯像元所得的遥感信息为Sveg，混合像元的植被成分所贡献的信息Sv可以表示为Sveg与FC的乘积。

同理，设全由土壤所覆盖的纯像元所得的遥感信息为Ssoil，混合像元的土壤成分所贡献的信息SS可以表示为Ssoil与1－FC乘积

将(4)、(5)两式代入(3)可得

其中Ssoil和Sveg都是参数，因而可以用上面式子计算植被覆盖度。

图1 密云县归一化植被指数植被图

图2 彩色合成方法与各颜色相应的解释

Bradley将归一化植被指数与植被覆盖度作线性相关分析.肯定了归一化植被指数与植被覆盖度有良好的相关性。此项研究表明归一化植被指数符合像元二分模型的条件，可以代入模型[9]。

由上述公式变换可得下面的利用归一化植被指数计算植被盖度的公式:

其中:NDVIsoil为完全裸土或无植被覆盖区域的值，即无植被像元的NDVI值;NDVIveg则代表全植覆盖像元的归一化植被指数值，即纯像元的归一化植被指数值。

根据公式 FC=(NDVI－NDVIsoi)/(NDVIveg－NDVIsoil)计算植被盖度，这里面最重要的是关于参数NDVIsoi和NDVIveg的选择。理论上来说NDVIsoi和NDVIveg应该为0和1，但在实际的应用当中，由于大气影响及地表湿度、粗糙度、土壤类型、土壤的影响，NDVIsoi会随着空间而变化，变化范围一般在－0.1～0.2之间[10];本文采取置信度为5%，在这个区间内取最大值和最小值。1992年、2005年及2009年的最大值和最小值分别为(0.588，0.053)，(0.508，－0.02)，(0.625，0.004)。利用 ENVI进行计算得到植被盖度等级图(图3)。

图3 植被盖度等级图

3.4 NDVI合成图及植被盖度等级合成图

将1992、2005、2009年的归一化植被指数图及植被盖度等级图分别进行(RGB)假彩色合成。根据RGB假彩色合成原理可知，当3幅影像的像元值都比较大时，像元呈灰白色;当像元值比较小时，像元呈黑色;其余颜色像元值变化比较大，为植被不稳定区域(图4和图5)。

图4 NDVI合成图

图5 植被盖度等级合成图

图6 NDVI监督分类图

3.5 监督分类

通过实地的抽样调查，配合人工目视判读，对遥感图像上抽样区的影像地物类别属性拥有先验知识，进行训练样区的选择，然后利用ENVI遥感图像处理软件中最大似然法对归一化植被指数合成图和植被盖度等级合成图进行监督分类。

图7 植被盖度等级监督分类图

依据密云县土地利用现状图，森林资源一类清查资料及Google Earth将植被盖度划分为3个等级，低植被盖度(0≤fc＜0.3)、中植被盖度(0.3≤fc＜0.7)和高植被盖度(0.7≤fc＜1)3 种类型。将三期图像进行RGB假彩色合成，(ABC)、(abc)和(123)分别代表1992、2005和2009年低、中和高植被盖度。结合假彩色合成原理，将两幅监督分类图进行统计，结果如表1和表2所示。

分类结果分析。参考实地考察密云县二类调查表，采用随机抽样的方法在研究区抽样，每期取点50个左右，利用相关的遥感软件计算Kappa指数。结果表明，在图6的8个类别当中，总体分类精度达85.25%，黑色 Kappa指数为1，品红色 Kappa指数最低为0.6175，其余均在0.7 ～0.8。在图7 的8 个类别当中，总体分类精度达83.13%，黑色Kappa指数为1，黄色Kappa指数最低为0.6532，其余均在0.68 ～0.85。

4 结论与讨论

2009年时密云县水域面积为6671.92 hm2，占全县面积比例是3.02%;密云县建筑用地面积4067.01 hm2，占全县面积比例是 1.82%;品红色区域代表的是从1992年到2005年这段时间内，归一化植被指数下降，从2005年到2009年这段时间，归一化植被指数升高，这部分区域主要是封山育林区域和农田植被以及部分建筑用地。面积为66806.37 hm2，所占比例为29.82%;黄色区域是2005年以后植被指数下降区域，主要是一些新开垦的荒地和裸地，及植被受到破坏的未利用地，面积为95543.82 hm2，所占比例为42.64%;和1992年相比2009年的时水域面积减少9947.9424 hm2，占全县面积比例的4.44%，这部分在图中的颜色是青色、蓝色和绿色;灰白色为植被稳定的区域，这部分区域面积40902.21，所占密云县面积的比例为18.26%。

表1 NDVI合成图监督分类结果

表2 植被盖度合成图监督分类结果

像元分类后对比法中各颜色所代表的意思和用逐个像元分类对比法略有不同，黑色代表的是低植被盖度区域，主要是水域和居民地，这部分面积 11584.71 hm2，所占比例为 5.29%;红色区域代表的随着时间的推移是该区域从高植被盖度向低植被盖度转移，主要是在这17 a中新增加的建筑用地和工矿用地。品红色和黄色所代表的意思和逐个像元分类对比法中的一致，都是植被不稳定的区域。逐个像元分类对比法中，像元植被不稳定区域所占比例为72.46%，在植被盖度分类后对比中所占比例为70.67%，可以看出两种方法中，所得出植被不稳定的区域相差不大。青色区域是无植被或者低植被覆盖区域向高植被覆盖区域转化，1992到2005蓝色区从低植被向中植被盖度发展，从2005年以发展到高植被盖度的。绿色区域是1990年到2005年从低植被盖度向高植被盖度发展，从2005年以后植被盖度下降到中植被盖度，这3种颜色主要是密云水库水域面积减少而引起的。相比分类后水域面积减少9634.26825 hm2，占4.3%，归一化植被指数逐个像元对比法更准确。

结合实际的调查和相关资料，密云县的这种植被动态变化与实际情况相符合，本文中所采用的两种方法都可以用于植被的大尺度监测，具有简便，灵活，操作性强等特点，使大范围内植被的监测变得简便易行而且准确。但逐个像元分类对比的方法只能用于3个分量的归一化植被指数彩色合成，所以对于长时间序列数据的分析处理，这是该方法的局限性。

[1]周文佐，潘剑君.南京市城市绿地现状遥感分析[J].遥感技术与应用，2002，17(1):22－26.

[2]章文波，符素华，刘宝元.目估法测量植被覆盖度的精度分析[J].北京师范大学学报:自然科学版，2001，37(3):402－408.

[3]张云霞，李晓兵，陈云浩.草地植被覆盖度的多尺度遥感与实地测量方法综述[J].地球科学进展，2003，18(1):85－93.

[4]赵春玲，李志刚，吕海军，等.中德合作宁夏贺兰山封山育林育草项目区植被覆盖度监测[J].宁夏农林科技，2000(6):6－14.

[5]Zhou Q，Robson M，Pilesjo P.On the ground estimation of vegetation cover in Australian rangelands[J].International Journal of Remote sensing，1998，19(9):1815－1820.

[6]张光辉，梁一民.黄土丘陵人工草地径流起始时间研究[J].水土保持学报，1995，9(3):78－83.

[7]黎夏.利用遥感与GIS对农田损失的监测及定量评价方法:以东莞市为例[J].地理学报，1997，52(3):279－287.

[8]李苗苗，吴炳方，颜长珍，等.密云水库上游植被覆盖度的遥感估算[J].资源科学，2004，26(4):153－159.

[9]Bradley C Rundquint.The influence of canopy green vegetation fraction on spectral measurements over native tall grass prairie[J].Remote Sensing of Environment，2002，81(1):129－135.

[10]Ro Hin A，Maspimave J，Lafleur J.Pressure dro through nonwoven geotex tiles:a new analytical model[J].IEEE International Conferenee on Geotextiles，1982，3:161－166.

[11]李本纲，陶澍.一种利用多时相TM影像分析地表植被的新方法:以敦煌地区绿洲植被变化为例[J].遥感学报，2000，11(4):295－298.

[12]刘广峰，吴波，范文义，等.基于像元二分模型的沙漠化地区植被覆盖度提取:以毛乌素沙地为例[J].水土保持研究，2007，14(2):268－271.

[13]陈巧，陈永福.Quick Bird遥感数据监测植被覆盖度的研究[J].林业科学研究，2005，18(4):375－380.

[14]Carlson T N，Ripley D A.On the relation between NDVI，fractional vegetation cover，and leaf area index[J].Remote Sensing of Enviroment，1997，62(3):241－252.

[15]Kaufman Y，Tanre D.Atmospherically resistant vegetation index(ARVI)for EOS-MODIS[J].IEEE Transactions on Geoscience and Remote Sensing，1992，30(2):261－270.