基于像元二分模型的植被覆盖度反演-以北京为例

1、基于像元二分模型的植被覆盖度反演-以北京市为例王玲（西北大学 城市与环境学院，陕西 西安710127）摘要：采用遥感技术监测植被覆盖度具有重要意义。本文以北京市为例，基于2013年的Landsat8 OLI影像，选取NDVI值为参数，采用像元二分模型对植被覆盖度进行反演，最终反演的结果与实际情况符合，说明采用该方法反演植被覆盖度可行。关键词：植被覆盖度、像元二分模型、NDVI、植被指数引言植被覆盖度（Vegetation fractional cover,简称fc）是指植被（包括叶、茎、枝）在地面 的垂直投影面积占统计区总面积的百分比，即植土比。通常林冠称郁闭度，灌草等植被称覆盖度1。它是衡量

2、地表植被覆盖的一个最重要的指标，被覆盖度及其变化是区域生态系统环境变化的重要指示，对水文、生态、全球变化等都具有重要意义2。根据监测手段，测量植被覆盖度的方法可分为地面测量和遥感测量两大类，测量常用于田间尺度， 遥感估算常用于区域尺度。目前已经发展了很多利用遥感测量植被覆盖度的方法，较为实用的地面测量单的方法就是目估法，缺点主要是主观性太强。客观的测量方法有样点法、样方法、样带法等， 借助于采样仪器的测量方法，空间定量计、移动光量计等。这些方法虽然提高了测量精度， 但野外操作不便，并且成本较高，难以在大范围内快速提取植被覆盖度。而采用遥感技术为监测大面积区域的植被覆盖度，甚至全球的植被覆盖度监

3、测提供了可能3。目前已有许多利用遥感技术测量植被覆盖度的方法，其中应用最广泛的方法是利用植被指数近似估算植被覆盖度，常用的植被指数为 NDVI。 一、数据源本文选取两景覆盖北京市的Landsat8 OLI影像、土地覆盖类型图以及北京行政边界矢量数据为数据源。其中，土地覆盖类型图是作为掩膜文件使用，其目的是为了便于植被覆盖度的估算；北京行政边界矢量数据是为了将两景镶嵌好的影像数据进行裁剪使用的，其目的是裁剪出北京市行政区内的范围。另外，Landsat8 OLI影像是从地理空间数据云网站上下载得到的，其成像时间为 2013年10月份。与Landsat7的ETM+成像仪相比，OLI成像仪获取的遥 感

4、图像辐射分辨率达到12比特，图像的几何精度和数据的信噪比也更高。OLI成像仪包括9个短波谱段（波段1波段9）,幅宽185km,其中全色波段地面分辨率为15m,其他谱段地面分辨率为30m4。口-您地后情力ETM*度段名善中前St*KRCttMhum(mJ1ml日国mF Ccsstal308aM 2事0450451530口 闻 1 Blu*和(152D 周030Sand 2 Green0 52M),fi0530Band 4 Red。后330侬0303a nd 3 RedO'.G3IM.69O和的中甯4 py|R0 77MJW0301 560M65Q30Eland 53WE1 55D-1./

5、50期2 f00-2.3003C&同 7 SWR 22.Q9&-2.35030aMSPan15Bana a. PanO WDHS 90015Band NCittue加1表1 Landsat8 OLI陆地成像仪波段参数二、研究方法本文反演植被覆盖度所采用的是像元二分模型方法，像元二分模型是一种简单实用的遥感估算模型，它假设一个像元的地表由有植被覆盖部分地表(Sv)与无植被覆盖部分地表(S)组成，而遥感传感器观测到的光谱信息(S也由这2个组分因子线性加权合成，各因子的权重是各自的面积在像元中所占的比率，如其中植被覆盖度可以看作是植被的权重。因此，像元二分模型的原理如下： 遥感传感器

6、观测到的光谱信息(S)由有植被覆盖部分地表(Sv)与无植被覆盖部分地表(S)组成，可得出：S = S + SS公式 1假设一个像元中有植被覆盖的面积比例为 fc ,即该像元的植被覆盖度，则裸土覆盖的面积 比例为1 -fc ,如果全由植被所覆盖的纯像元所得的遥感信息为Sveg ,则混合像元的植被部分所贡献的信息SM以表示为Sve冉fc的乘积：Sv =fc Sveg公式 2那么，Ss =(1 -fc )- Ssoil公式 3将公式2与公式3代入到公式1中，可得到：S =fc - Sveg +(1 -fc)S soil 公式 4对公式4进行变换，可得以下计算植被覆盖度的公式：fc =(S -Ssoi

7、l) (Sveg -Ssoil ) 公式 5其中Ssoil为纯土壤像元的信息，Sveg为纯植被像元的信息，因而可以根据公式5利用遥 感信息来估算植被覆盖度。将归一化植被指数(NDVI)代入公式5可以被近似为：fc =(NDVI - NDVIsoil ) (NDVIveg -NDVIsoil) 公式 6其中，NDVIsoil为裸土或无植被覆盖区域的NDVI值，即无植被像元的NDVI值;而2 NDVIveg则代表完全被植被所覆盖的像元的NDVI值，即纯植被像元的NDVI值。当区域内可以近似取 VFCmax=100%, VFCmin=0% VFC= (NDVI -NDVImin)/ ( NDVIma

8、x -NDVImin), NDVImax和NDVImin分别为区域内最大和最小的NDVI值。由于不可避免存在噪声，NDVImax和NDVImin一般取一定置信度范围内的最大值与最小值，置信度的取值主要根据图像实际情况来定；当区域内不能近似取VFCmax=100% VFCmin=0%当有实测数据的情况下，取实测数据中的植被覆盖度的最大值和最小值作为VFCma坏口 VFCmin,这两个实测数据对应图像的NDVI作为NDVImax和NDVImin。当没有实测数据的情况下，取一定置信度范围内的NDVImax和NDVImin。VFCma/口VFCmin艮据经验估算。 三、数据处理 1、数据预处理本文使用

9、的Landsat8 OLI为L1T级别数据，不需做几何校正处理。而北京市需要两景 Landsat OLI数据覆盖，因此首先要进行图像镶嵌和裁剪，然后进行大气校正等预处理过程。(1) 辐射定标辐射定标是将传感器记录的电压或数字值转换成绝对辐射亮度的过程。目的是消除传感器本身所产生的误差，由于传感器在不断的运行中光学器件性能逐渐退化，因此定标的系数也随之不同，这些定标系数也在不断的更改，在用户获得数据的时候，这些定标系数也在影像的头文件中同时提供给用户。遥感数据辐射定标就是将传感器得到的灰度值转换为星上的辐射亮度值或星上反射率，即表观辐射度或表观反射率。辐射定标主要校正由传感器的灵敏度 带来的辐射

10、误差 。其目的是为FLAAS次气校正准备数据：定标符合单位要求的辐射量数据、转换数据储存顺序等。该处理过程在 Envi5.2中实现，具体操作：在 ENVIToolbox中，选择Toolbox/Radiometric Correction/ Radiometric Calibration ,选择*_MultiSpectral 多光谱组(7个波段),打开辐射定 标工具，对两景影像分别做辐射定标。(2)影像镶嵌因本文所使用的影像数据源是两景Landsat OLI影像，因此需进行影像镶嵌，镶嵌的目的是将不同的影像文件无缝地拼接成一幅完整的包含研究区域的影像。该处理过程在Envi5.2中实现，具体操作：

11、在Toolbox中，选择/Mosaicking/Seamless Mosaic ,打开无缝镶嵌 工具，然后进行相关参数设置，如下所示：(3)影像裁剪因本文所使用的影像数据包含了北京市行政区划以外的部分地区，因此需进行影像裁剪,以将研究区裁剪出来，并且减小了数据量，加快了数据处理速度，本文使用北京行政边界矢 量裁剪图像。过程在 Envi5.2中的具体操作如下：在Toolbox中，选择 /Regions of Interest/Subset Data from ROIs ,打开裁剪工具：影像裁剪结果如下所示:(4) Flaash大气校正，只电磁波在大气中的传输和遥感器观测过程中受光照条件以及大气作

12、用等的影响有小部分（在0.85um波段80%,在0.45um波段50%）太阳辐射能反射到遥感器,导致遥感器的测 量值与地物实际的光谱辐射率不一样。辐射损失主要发生在大气吸收和散射过程，因此地表参数的遥感定量反演研究中，必须纠正目标辐射的不确定性信息6。ENVI中的FLAASK型是基于MODTRAN4+射传输模型，通过参数查找表来进行大气校 正的商业化软件。FLAASH：气校正模块支持多种传感器数据，其光谱处理范围0.4 wm-2.5可以有效地去除水蒸气/气溶胶散射效应，同时该方法基于图像像素级的校正，能够解决目 标像元和邻近象元的“邻近效应”问题，校正结果精度高，简单易行7。然后，对大气校正前

13、后同一地物的光谱曲线进行对比，这里以植被为例，光谱曲线如下 图所示：校正前的植被光谱曲线校正后的植被光谱曲线2、植被覆盖度估算(1)计算 NDVI本文选取NDVI值为参数，采用像元二分模型对植被覆盖度进行反演，根据植被覆盖度的计算公式可知，要求取植被覆盖度，首先需要计算NDVI。在Envi5.2中的具体操作如下:在Toolbox中，选择 Spectral/Vegetation/NDVI , NDVI Calculation Input File 面板中，选择LC8_rad_beijing_ref.dat 图像，求算 NDVI,如下：XNDVI求算结果如下:由于大气校正后的结果有部分像元为负值，

14、主要集中在阴影地区，这部分区域计算得到的NDVI在卜1,1之外，为了便于后面的分析， 我们这里统一将这部分像元进行处理， 即NDVI 值大于1的变为1,小于-1的变成-1。在Bandmath中的表达式为：-1>b1<1 ,其中bl代 表的是NDV|得到去除异常值文件：NDVI上除异常值.dat(2) 掩膜文件制作该过程主要是为了计算NDVI的最大值、最小值所服务的，根据土地利用分类图(共5类，林地、农业用地、城市用地、水体与其他)制作各种土地利用类型的掩膜文件，在Envi5.2中的具体操作如下：在Toolbox中选择/Raster Management/Masking/Apply

15、Mask ,打开制作掩膜工具：采用该方法制作林地、农业用地、城市用地、水体与其他的掩膜文件，其中林地与耕地的掩膜文件制作结果如下:（2）获取阈值这一步就是求解NDVImax和NDVImin,使用上一步获取的掩膜文件分别对 NDVI图像文件 进行统计，在一定置信度范围内获取每个掩膜文件（也就是土地覆盖类型） 对应的最大和最小NDVI值。在Toolbox中，选择/Statistics/Compute Statistics ,进行统计，然后在统计结果中，取一定的置信度获取最大和最小的 NDVI值。如这里的林地覆盖区域的统计结果（如下图），这个过程带有很大的主观性，我们需要根据统计学原理自己制定一套规

16、则（比如5%的置信度），这里我就以NDVI值对应像元数量增加到 5位数字为置信区间，选择 NDVImin=0.3804 ,NDVImax=0.8667。0.356982361449590.10121.78120.364T06321535910.10611.66730.372592631628S40,11332.00110.3004 _10000LT26S40.12292.1240a. 3S3210879183T330.133T2. 25TT0.395L11579L953120.14232 40000. 4039121412074530.14922.54920.8431veneososoi0.

17、588T0.85103166880921000.309499 4362口.05882019181】2沏0.2431S9.68430.86671222081246000.150269 33450 3T456992S131592u.005959 9204C.03243&55SL35257. 0450取 9654同样的方法得到其他地物覆盖类型的NDVI阈值，其中，水体没有植被(水藻不属于植被)认为这部分区域的植被覆盖度为0,如下表：士地浮米类型NDVImin (NDVIsoil)NDVImax (NDVIveg)林地O.Kbtj；农业用地0.16040.7794班巾.用地0.04040.47

18、89水体ii：其他(3)生成参数文件植被覆盖度的计算公式：fc =(NDVI - NDVIsoil ) (NDVIveg -NDVIsoil)该过程是根据上面得 到的NDVI阈值分别生成NDVIsoil和NDVIveg参数文件，也即 NDVImin与NDVImax。该过程主要 使用Envi5.2的bandmath工具，并且：NDVIsoil: b1*0.3804+b2*0.1604+b3*0.0404+b4*0+b5*0.0946其中，b1:林地掩膜文件b2：农业用地掩膜文件b3:城市用地掩膜文件b4：水体掩膜文件b5:其他用地掩膜文件NDVIveg: b1*0.8667+b2*0.7794+

19、b3*0.4789+b4*0+b5*0.5347其中，b1:林地掩膜文件b2：农业用地掩膜文件b3:城市用地掩膜文件b4：水体掩膜文件b5:其他用地掩膜文件最终，生成的参数文件如下所示:(4)植被覆盖度估算利用上一步得到的 NDVIsoil和NDVIveg参数文件带入公式：fc =(NDVI - NDVIsoil ) (NDVIveg-NDVIsoil）,该过程也是利用 Envi5.2中的Bandmath工具来实现的，其表达式为：（b1-b2）/（b3-b2）, 其中，bl为NDVI（对应的文件名为 NDVI_去除异常值.dat"）、B2为NDVIsoil参数文件、B3:为 NDVI

20、veg参数文件，植被覆盖度估算结果如下：刁 qroprirks View He pLit： +_* r lt.an Lh： iib* w uiWj £市加生：j jiwth.偿-mB H 3- Yl riPmvisw HhJ ¥alh (IL'J>| Hnm*5 Jiti. 0 '.B 超在 12 0：-:U Bstd *U& C:1* J我们分析下结果，会发现有一些异常值，即值在0,1之外，这些异常值是在NDV置信度 之外的那部分像元产生的（也包括NDVI异常像元）。这些像元数量不多，大约占3.7%左右。还有背景和水体区域的植被覆盖度的值为-

21、NaN,即无效值，因为分母为 。造成的。第一种异常值可以将小于 。的值变成0,大于1的值变成1,用bandmath工具即可， Bandmath表达式为：0.0>b1<1.0，其中b1为植被覆盖度 ；-NaN可以用掩膜进行处理，即在Build Mask中用-NaN生成掩膜。去掉异常值之后，并对其进行分类显示，最终得到的植被覆盖度图如下：(5)结果验证目前业内植被覆盖度的验证方法主要是通过野外调查结果进行验证，如下为一种方法：以与地面垂直的角度用数码相机拍摄采样点的地面照片，使用GPS1位获得采样点的经纬度坐标。为获得准确的植被覆盖度实测数据减少像片边缘变形误差，将数码相机得到的数字图像截取长、宽各三分之二的中心地带，采用非监督分类法为10类，并将分类结果分为植被、非植被两类，以此来计算出植