基于遥感影像的北京市植被覆盖率反演.docx

基于遥感影像的北京市植被覆盖率反演.docx

《基于遥感影像的北京市植被覆盖率反演.docx》由会员分享，可在线阅读，更多相关《基于遥感影像的北京市植被覆盖率反演.docx（26页珍藏版）》请在冰豆网上搜索。

基于遥感影像的北京市植被覆盖率反演

基于ETM遥感影像的北京市植被覆盖度反演

1绪论

1.1研究区概况

北京中心位于北纬39度54分20秒，东经116度25分29秒。

位于华北平原西北边缘。

毗邻渤海湾，上靠辽东半岛，下临山东半岛。

北京与天津相邻，并与天津一起被河北省环绕。

西部是太行山山脉余脉的西山，北部是燕山山脉的军都山，两山在南口关沟相交，形成一个向东南展开的半圆形大山弯，人们称之为“北京弯”，它所围绕的小平原即为北京小平原。

诚如古人所言：

“幽州之地，左环沧海，右拥太行，北枕居庸，南襟河济，诚天府之国”。

全市平均海拔43.5米。

北京平原的海拔高度在20～60米，山地一般海拔1,000～1,500米。

北京市国土面积16410.54平方公里，市区面积12187平方公里，建成区面积1386平方公里。

1.2数据的选择

本研究的研究的是北京市植被覆盖度反演，应选择植被生长的旺季，考虑到北京市的气候和季节，应选择6~10月份之间的数据影响为最合适的。

据此，本研究选择的是2009年9月22日的landsat5的TM遥感影像。

图1-1图像的示例图：

影像绿地信息明显，有较好的识别效果。

基本上满足研究的需要。

1-1影像数据

1.3研究的意义

植被，包括森林、灌丛、草地和农作物，既是生态系统的主要组成部分，也是生态系统存在的基础，具有截流降雨、减缓径流、防沙治沙、保持水土等功能，联结着土壤、大气和水分等自然过程，在陆地表面的能量交换、生物地球化学循环和水文循环等过程中扮演着重要角色，是全球变化研究中的“指示器”[1]。

植被根据生态系统中水、气等的状况，调控其内部与外部的物质、能量交换。

植被覆盖与气候因子关系极为密切，研究植被覆盖变化 对气候的影响是气候变化研究的主要内容之一，它影响着土壤湿度、地表温度和地表能量与水的循环。

一个城市的植被覆盖度可以反映一个城市规划的情况，间接的反应一个城市的环境质量，是城市的重要组成部分。

城市化的迅速推进，带来了多样化的生态足迹，植被覆盖度，土壤污染率，地表侵蚀率，逐渐成为生态研究的热点，也成为环境保护的重点。

2基本原理及反演技术流程

2.1植被覆盖度

植被覆盖度是指植被（包括叶、茎、枝）在地面的垂直投影面积占统计区总面积的百分比。

容易与植被覆盖度混淆的概念是植被盖度，植被盖度是指植被冠层或叶面在地面的垂直投影面积占植被区总面积的比例。

两个概念主要区别就是分母不一样。

植被覆盖度常用于植被变化、生态环境研究、水土保持、气候等方面。

植被覆盖度的测量可分为地面测量和遥感估算两种方法。

地面测量常用于田间尺度，遥感估算常用于区域尺度。

目前已经发展了很多利用遥感测量植被覆盖度的方法，较为实用的方法是利用植被指数近似估算植被覆盖度，常用的植被指数为NDVI。

2.2象元二分法

像元二分模型是一种简单实用的遥感估算模型，它假设一个像元的地表由有植被覆盖部分地表与无植被覆盖部分地表组成，而遥感传感器观测到的光谱信息也由这2个组分因子线性加权合成，各因子的权重是各自的面积在像元中所占的比率，如其中植被覆盖度可以看作是植被的权重。

下面简单介绍下在像元二分模型的基础上研究的模型：

………………………………………………（2-1）

其中,NDVIsoil 为完全是裸土或无植被覆盖区域的NDVI值，NDVIveg 则代表完全被植被所覆盖的像元的NDVI值，即纯植被像元的NDVI值。

两个值的计算公式为：

………………………（2-3）

……………（2-4）

利用这个模型计算植被覆盖度的关键是计算NDVIsoil和NDVIveg。

这里有两种假设：

1）当区域内可以近似取VFCmax=100%，VFCmin=0%。

公式（2-1）可变为：

………………………………………………（2-4）

NDVImax和NDVImin分别为区域内最大和最小的NDVI值。

由于不可避免存在噪声，NDVImax和NDVImin一般取一定置信度范围内的最大值与最小值，置信度的取值主要根据图像实际情况来定。

2）当区域内不能近似取VFCmax=100%，VFCmin=0%

当有实测数据的情况下，取实测数据中的植被覆盖度的最大值和最小值作为VFCmax和 VFCmin，这两个实测数据对应图像的NDVI作为NDVImax 和NDVImin。

2.3植被指数NDVI

本研究选择NDVI估算植被覆盖度，尽管NDVI增强了近红外与红色通道反射率的对比度，增强了低值部分，导致对少植被区域较高、对多植被取区域较低的敏感度，但有其自身的优点：

1）NDVI是植被生长状况的最佳因子。

2）经过比值处理的NDVI可以消除与太阳高度角、地形、云等和大气条件有关的辐射条件变化的影响，使在研究区植被覆盖度提取中的误差一定程度上地减小。

3）北京市的几种典型的地面覆盖类型在NDVI图像上区分明显，植被突出。

鉴于以上的优点，本研究用NDVI提取北京市的植被覆盖度。

NDVI，即近红外波段NIR与可见光红波段RED数值之差与这两个波段数值之和的比值即：

NDVI=

……………………………………………（2-5）

式中：

NIR为近红外波段数值；RED为可见光波段数值。

2.4主要技术流程

数据的获取

数据的预处理

辐射定标

大气校正

影像拼接

图像裁剪

监督分类

NDVI

NDVI

异常值处理

林地

农用地

城镇

水体

其他

掩膜文件

统计最大最小值

Lmin=b1*lmin+b2*lmin+b3*lmin+b4*lmin+b5*lmin

Lmax=b1*lmax+b2*lmax+b3*lmax+b4*lmax+b5*lmax

VFC=（b1-b2）/（b3-b2）

处理异常值后的NDVI

植被覆盖度

异常值去除

分级显示

图2-1技术路线

3数据的预处理

3.1数据的打开

本研究中图像为美国landsat5卫星所拍摄的TM5数据，包含有7个波段，格式为TIF。

带有地理坐标。

因本研究所使用的软件为ENVI5.0,在其中可以直接打开，跳过了以往的需要进行波段的叠合过程，具体过程为：

OpenAs→Landsat→GeoTIFFWithMetadata打开影像如下：

图3-1图像（01）

图3-2图像（02）

3.2影像的辐射定标

随着遥感技术的迅速发展,遥感应用逐渐从定性走向定量,辐射定标是实现遥感信息定量化的基本前提。

辐射定标是指建立遥感传感器的数字量化输出值DN与其所对应视场中辐射亮度值之间的定量关系。

通常遥感数字图像给出的是像元DN值。

利用DN值，只能进行同景图像内部的相对比较。

全球资源和环境变化研究要求遥感技术能够提供长时相、多区域、多种传感器的遥感数据。

只有将图像DN值转换成对应像元的辐射亮度值,才能对不同地点、不同时间和不同类型传感器获取的遥感数据进行定量比较与应用，以满足全球资源和环境变化研究的需要，而这个转换过程就称为辐射定标。

目前辐射定标包括3个方面的内容:

1）发射前的实验室定标;

2）基于星载定标器的飞行中定标;

3）在轨运行期间采用基于陆地（或海面）特性的“替代定标”,或借助其他卫星进行的“交叉定标”。

辐射定标贯穿卫星整个生命周期,是保持卫星数据精度的一项非常重要的基础工作。

通过辐射定标，可以监测传感器的性能变化，并定期或不定期地给出传感器的辐射定标系数。

定量遥感技术的发展，全球资源和环境变化的遥感监测以及多光谱、多时相和多种卫星传感器遥感数据的综合应用和定量分析技术的发展，迫切地对卫星传感器的辐射定标提出高精度的要求。

对于landsat5的TM影像的辐射定标是对于像元DN值的线性函数的变换，可以通过3-1公式进行模型的变换。

…………………………………………………（3-1）

式中：

L指定标后的反射率值；单位（mW/（cm2*Sr*

））

gain指增益系数；

B指偏移数值；

本研究是基于软件ENVI5.0，选择定标辐射工具后，选择相应的影像数据，弹出如下对话框：

3-3辐射定标

其中定标类型选择Radiance，输入的参数类型为BIL,数据的类型必须为Float,然后点击对话框中的FLAASHSettings按钮。

选择数据的输入路径。

3.3影像FLAASH大气校正

3.3.1大气组成及FLAASH大气矫正

大气主要包括气体分子和一些其他微粒。

分子主要有N2和O2，约占99%，其余的1%是03、CO2、H2O及其他的气体。

大气的微粒主要是烟、尘埃、雾霾、小水滴及气溶胶。

气溶胶是一种固体和液体的悬浮物，一般直径大小在0.01~30

之间，多分布在高度在5KM以下。

大气层自下而上依次分为对流层、平流层、中间层、热层和散逸层。

遥感利用的一切辐射都必须经过地球的大气层。

由于传感器的工作波段选择的原因，对于传感器接收到的电磁波辐射影响最大的是对流层和平流层。

当太阳辐射穿过大气时，必然要受到大气散射、吸收和折射的影响。

大气校正的目的是消除大气和光照等因素对地物反射的影响，获得地物反射率和辐射率、地表温度等真实物理模型参数，用来消除大气中水蒸气、氧气、二氧化碳、甲烷和臭氧对地物反射的影响，消除大气分子和气溶胶散射的影响。

FLAASH 可以处理任何高光谱数据、卫星数据和航空数据（860nm/1135nm），这些数据是由HyMAP、AVIRIS、CASI、 HYDICE、HYPERION（EO-1）AISA、HARP、DAIS、Probe-1、TRWIS-3、SINDRI、MIVIS、OrbView-4、NEMO等传感器获得的。

FLAASH还可以校正垂直成像数据和侧视成像数据。

Flaash大气校正使用了 MODTRAN4+ 辐射传输模型的代码，基于像元级的校正，校正由于漫反射引起的连带效应，包含卷云和不透明云层的分类图，可调整由于人为抑止而导致的波谱平滑。

FLAASH可对Landsat,SPOT,AVHRR,ASTER,MODIS,MERIS,AATSR,IRS等多光谱、高光谱数据、航空影像及自定义格式的高光谱影像进行快速大气校正分析。

能有效消除大气和光照等因素对地物反射的影响，获得地物较为准确的反射率和辐射率、地表温度等真实物理模型参数。

FLAASH大气校正对图像有以下几个要求：

1）数据是经过定标后的辐射亮度（辐射率）数据，单位是：

；

2）数据带有中心波长值，如果是高光谱还必须有波段宽度，这两个参数都可以通过编辑头文件信息输入；

3）数据类型----支持四种数据类型：

浮点型、长整型、整型和无符号整型；

4）数据的储存类型是ENVI标准栅格格式文件，且必须是BIL或者BIP;

5）波谱范围为：

400——2500nm

3.3.2大气校正过程

大气校正的影响必须是经过辐射定标后的影响数据，首先打开辐射定标后的影像数据，选择FLAASH大气校正命令，弹出如图3-4的对话框：

3-4图FLAASH大气校正

点击InputRadianceImage时会弹出如3-5的对话框，此对话框的目的是用来做单位转化的，因为在上面的辐射定标的时候已经进行过单位的转化了，再次不需要再次转化了，直接点击OK。

图3-5单位转化

FLAASH大气校正需要输入较多的参数，如影像的中心坐标、传感器类型、能见度、影像拍摄时间等参数。

这些参数均可以在头文件里查询。

输入各种参数之后以及储存路径之后点击OK。

3.3.3FLAASH大气校正前后比较

FLAASH大气校正之后基本上消除了大气散射等因素的影响，主要的差别是影像在光谱曲线图上的直观表达，通过查看相同区域的光谱曲线以及根据地物的光谱曲线可以检验FLAASH大气校正的结果的效果。

如图3-6和3-7是植被光谱曲线FLAASH大气校正前后的光谱曲线的对比。

图3-6大气校正前的植被光谱曲线图

图3-7大气校正后的植被光谱曲线图

图3-8植被的光谱曲线图

从上两光谱曲线图，结合图3-8植被地物光谱曲线图可以发现大气对于影像的影响经过FLAASH大气校正之后有了显著的改善，能够基本上和植被的光谱曲线吻合。

影像的质量上有了很大的提高，为下一步的植被覆盖度反演奠定了反射率上的基础。

3.4影像的拼接

由于北京市一幅影像的不能覆盖完全，需要两幅影像进行拼接，在两幅影响进行大气校正之后需要对两幅影像进行拼接。

因为landsat5TM影像有自己的地理坐标，所以我们在进行拼接的时候选择基于地理坐标的拼接，其过程如下：

打开图像拼接命令，选择经过大气校正后的两幅影像，如图3-9。

图3-9影像的拼接

在上层的影像右击，选择EditEntry命令，弹出如图3-10对话框，在DataValuetoIgnore中输入0，点击OK。

图3-10EditEntry命令对话框

点击file菜单，选择Apply，弹出如下对话框（图3-11），选用三次卷积从采样方式，选择输入路径，背景值输为0，点击OK。

即完成图像的拼接。

图3-11图像的拼接对话框

拼接的结果如图3-11.

图3-12影像拼接后

3.5影像的裁剪

本研究的区域为北京市，两幅影像不仅包含北京市而且还有北京市边缘地区，属于非研究区域，在进行植被覆盖度反演的时候会影像数据的处理，而且会造成数据的冗余。

因此，将研究区以外的区域裁剪掉对于研究的精度和效率在一定程度上有所提高。

本研究是基于软件ENVI5.0，改进了以往需要将边界矢量数据转换成栅格数据，因此，可以在拼接的基础上直接打开矢量数据边界，其裁剪的技术路线如图3-13

影像图

裁剪后影像图

边界矢量数据

图3-13影像裁剪路线

裁剪的过程如下：

打开拼接后的影像，同时打开北京市的矢量边界文件，矢量边界会自动的叠加到影像上如图3-14.

图3-14影像边界叠加

选择裁剪命令，在弹出的对话框中选择适量边界文件，并将去背景值设置为0，选择输出路径，点击OK，完成影像的裁剪，如图3-15。

图3-15裁剪后的影像

4.影像的信息提取

4.1计算NDVI及异常值处理

4.1.1计算NDVI

本研究是利用NDVI值进行植被覆盖度的反演，因此要对影像进行计算NDVI值。

由公式2-5可知，利用影像的第三波段和第四波段进行求解。

计算过程如下：

打开NDVI命令，在弹出的对话框选择裁剪后的北京市影像，点击OK。

在弹出的对话框中选择输入的影像类型为landsatTM，在波段选择里选3和4波段并选择输入的路径，点击OK，完成NDVI的计算。

结果如图3-16.

图4-1NDVI

4.1.2NDVI异常值处理

影像因为受到太阳高度角、地形、阴影等因素的影响，很少一部分区域计算后的NDVI值会出现大于1和小于0的情况，属于异常值，在进行植被覆盖度反演之前要进行处理。

本研究中将小于0的NDVI值归并为背景值，即赋值为0，将大于1的赋值为1。

处理的方法是用波段运算，处理过程如下：

打开波段运算命令，输入计算表达式（b1lt0）+（b1gt1）,点击0K。

查看异常值的区域如图3-17

图4-2NDVI异常值

其中的一个的大区域的是水体，可以看出符合正常误差范围里。

打开波段运算命令，输入计算表达式：

（b1lt0）*0+（b1gt1）*0+（b1ge0andb1le1）*b1，其中b1定义为NDVI值，选择输入的路径，点击ok，完成异常值的处理，结果如下图3-18

图4-3异常处理后的NDVI

4.2用监督分类进行土地分类

本研究是利用像元二分法进行植被覆盖度的反演，其中主要的一个步骤就是利用监督分类方法对影像进行土地的分类，然后统计不同土地类型的NDVI值的最值进行植被覆盖度的反演，因此在反演之前应进行监督分类进行土地的分类。

研究中，把北京市的土地分成五大类，即林地、农用地、城市用地、水体、其他用地。

利用FLAASH大气校正后的影像进行监督分类，过程如下:

在ENVI5.0中，监督分类是一个工作流，打开分类工作流，选择裁剪后的影像，添加五类用地并选择训练样本，如图4-4

图4-4监督分类训练样本

最终的分类结果如图4-5，

图4-5土地分类图

4.3利用土地分类图及NDVI统计各类的最值

本研究用的方法是利用像元二分法进行植被的覆盖度反演，根据公式2-1可知需要统计每一类土地的NDVI最小值近似的代替NDVIsoil，统计出每一类土地NDVI的最小值代替植被纯像元的NDVIvge。

根据各类土地的最值求出植被覆盖度VFC。

因为每一类的都需要统计，如果对每一类土地都进行求解，不仅繁琐而且容易出错，所以利用波段运算。

计算公式如下：

…………………………….（4-1）

式中：

Lmin是影像所有土地类别的最小NDVI值（是面积上的统计，并非数值的求和）；

b1、b2、b3、b4、b5是土地的类别；

lmin是每一类土地的NDVI的最小值；

………………………….（4-2）

式中：

Lmax是影像所有土地类别的最大NDVI值（是面积上的统计，并非数值的求和）；

b1、b2、b3、b4、b5是土地的类别；

Lmax是每一类土地的NDVI的最大值；

4.3.1利用土地分类图的掩模制作

如果要统计每一类土地的NDVI的最值就需要讲每一类单独分类出来进行统计，因此需要通过制作掩膜文件。

本研究中是根据分类图中的每一类别的分类值进行掩模的制作，分类出每一土地类别的NDVI值进行统计。

过程如下:

1）打开经过异常值处理的DNVI影响图和分类图，打开ComputeStatistics命令，在弹出的对话框中选择经过异常值处理的NDVI影像，点击MsakOption下拉菜单，选择BuidlMask命令，在弹出的对话框中选择Option下来菜单；

2）因为需要利用分类图进行掩模的制作，而分类图对应的是每一类的数值，所以选择ImportImageRange，在弹出的对话框中选择分类图影像，点击OK，弹出图4-6对话框，

图4-6ImportDataRange

3）在4-6对话框中输入对应土地的数值，然后点击ok，完成一次掩膜的制作，图4-7为林地的掩膜文件。

图4-7林地的掩膜影像文件

4）重复上述的步骤，完成各类土地的掩膜文件。

4.3.2NDVI的最值统计

在完成上述的掩模文件的制作之后，就需要对各类的土地的NDVI的最值进行统计，为之后的植被覆盖度的反演做基础。

在本研究的的植被覆盖度反演公式中理论上是针对每一个像元的植被和裸地的NDVI值，但在实际中，因为不同时期植被长势不同以及受邻近的像元的影响，NDVI的最值要根据研究区的情况而定。

本研究取统计直方图的第一个异常点所对的NDVI值作为每一类的最值进行植被覆盖度的反演。

具体过程如下：

1）打开ComputeStatistics命令，在弹出的对话框中选择进行异常值处理后的NDVI影响图，点击SelectMaskBand，选择制作好的掩膜文件，点击OK。

2）在上述弹出的对话框中勾选直方图前的复选框，同时输出直方图和统计表格，如图4-8，

图4-8统计结果

3）重复上述的步骤，完成每一类的土地NDVI的最值统计，统计结果如表4.1。

表4.1DNVI最值统计

土地类别

统计最小值

统计最大值

林地

0.4862

0.9725

农用地

0.3298

0.8093

城镇用地

0.0848

0.5206

水体

0

0

其他用地

0.0227

0.5508

4.4植被覆盖度反演

4.4.1植被覆盖度反演计算

根据公式4-1和4-2以及植被覆盖反演公式，需要求出Lmin和Lmax来求解植被覆盖度，因此，首先要对影像的最值求解，然后进行反演。

具体的过程如下：

1）打开bandmath命令，在弹出的对话框中输如下公式：

b1*0.4862+b2*0.3298+b3*0.0848+b4*0+b5*0.0227点击ok。

2）在弹出如下的对话框中分别定义b1、b2、b3、b4、b5为林地、农用地、城市用地、水体、其他用地所生成的掩模文件，选择输入路径然后点击ok，完成最小值的生成。

图4-9最小值

3）重复上述步骤，在对话框中输入公式：

b1*0.9725+b2*0.8093+b3*0.5206+b4*0+b5*0.5508完成最大值的文件统计。

4）根据本研究中公式2-1计算植被覆盖度的反演，同时为了让背景值不参与运算，因此在进行植被覆盖度的反演时需要对公式进行变化。

5）打开bandmath命令，在弹出的对话框中输入公式：

（B2ne0.0）*（b1-b2）/（b3-b2）点击ok。

6）在打开的对话框中分别需要定义b1、b2、b3分别是去除异常值的NDVI、最小统计影像、最大影像值影像，选择输出的路径然后点击ok，完成植被覆盖度的反演。

7）在反演的结果中，因为有极少部分的像元会出现超出执行范围之外（理论上反演结果应该在0~1之间），因此需要对结果进行归一化处理。

8）同样的打开bandmath命令，在弹出的对话框中输入表达式：

0.0>b1<1.0点击ok，在弹出的对话框中定义b1为上述反演结果影像图，选择输出的路径，点击ok。

如图4-10，

图4-10植被覆盖度影像

4.4.2植被覆盖度反演影像分级显示

在完成植被覆盖度反演的结果中，为了能够直观并且美观的看出区域的植被覆盖度的情况。

通常对结果影像进行分级的显示，具体步骤如下：

1）在ENVI5.0的layermanger中选中最终反演结果，右击鼠标，选择影像分级命令，弹出如下图对话框，其中默认16级，本研究中我们分成十个级别，分别定义其区间和颜色，如下图，

图4-11影像分级显示

2）完成上述的定义区间和颜色的选取，点击ok，完成影像植被覆盖度的反演分级，如下图为分级结果，

图4-12分级显示影像图

5结论

本文对植被覆盖度反演的研究是基于像元的二分法，通过对北京市的区域的影响进行具体的分析与研究，对植被覆盖度反演有了一定程度的认识和理解。

此次研究翻阅大量的相关资料，深刻理解遥感技术在城市规划研究和城市土地利用等方面的优势，当然某些程度上的不足。

遥感技术在理论上存在不足，还需要不断的深入和探究。

至此也要感谢王召海老师的大力指导和帮助，在理论和技术上细心地的讲解和指导。

同时也感谢同学的热心的帮助。