基于MODIS数据的北京西北部土地资源监测研究

(1.国土资源部土地利用重点实验室，中国土地勘测规划院，北京，100035；2.中国测绘科学研究院，北京，100039)

本文以MODIS 16天NDVI时间序列数据、8天LST数据、1 ∶ 5万DEM数据等辅助数据为基础，主要论述了京西北地区土地资源调查和多年土地利用与植被覆盖变化研究。首先，选择适合MODIS数据分类的土地覆盖分类系统，利用主成分分析法对NDVI时间序列数据进行信息增强和压缩。结合LST数据、DEM数据和降雨温度数据，采用模糊K-means无监督分类方法对研究区进行土地覆盖分类，得到土地资源现状。然后，运用CVA分析方法，分析了京西北地区土地利用和植被覆盖的多年变化。结果表明，MODIS数据能很好地应用于大范围的土地资源监测，并能取得较好的效果。

关键词:北京西北；莫迪斯；土地资源现状；土地利用和植被变化

随着人口、资源和环境之间的矛盾日益尖锐，为了实现可持续发展的战略目标，世界各国政府都在大力加强资源和生态环境监测系统的建设。我国《全国生态环境建设规划》和《全国生态环境保护纲要》也明确提出要完善生态环境监测和信息服务体系。国土资源部“十五”科技发展规划纲要强调，大力推进国土资源管理信息化，努力实现国土资源调查评价现代化。其中，土地资源调查与监测是其主要内容之一。

随着现代遥感技术的飞速发展，各种适用于不同空间尺度的土地资源调查与监测的遥感数据相继出现。目前，SPOT5、Landsat TM/ETM ++遥感数据是土地资源调查的主要数据源，适用于乡、县、区等不同尺度，但应用于更大空间尺度的土地资源调查时，耗费大量人力物力，且不经济。近年来中空间分辨率和高时间分辨率的MODIS数据为大比例尺土地资源调查提供了较好的数据源。Tucker等人的研究表明[1]，归一化植被指数NDVI实际上反映了植被生物量、盖度和叶绿素含量的生物物理化学性质，不同相位条件下的NDVI序列可以准确反映植被生长的季节变化规律。这已经成为利用遥感数据进行大面积植被和土地覆盖制图的基本思路[2 ~ 4]。

本文试图以MODIS NDVI时间序列数据集为主要数据源，结合MODIS LST、DEM、降雨量和温度等辅助数据，首先选择合适的土地覆盖分类系统，并采用模糊K-means无监督分类方法对京西北地区的土地覆盖自动分类进行研究。然后利用变化向量(CVA)分析方法分析该区土地利用和植被覆盖的多年变化，为大尺度土地资源调查和监测提供一种快捷方便的方法。

1研究区概况

北京西北地区是中国生态环境建设部门关注和投资的重点。研究区主要包括北京西北部的风沙源区，涉及河北、山西、内蒙古8个市(州、盟)51个县(市、旗)，总土地面积22.83×104 km2。该地区西起内蒙古四子王旗，东至内蒙古敖汉旗，南至山西代县，北至内蒙古阿鲁科尔沁旗。地理坐标为东经110 20 ' ~ 121 01 '，北纬38 565438。

研究区位于内蒙古高原中部，黄土高原北端，内蒙古、山西、河北三省交界处。区内地表形态主要由高原、山地、丘陵和盆地组成，地势中间高，南北低。研究区横跨中温带和寒温带，属于干旱半干旱大陆性季风气候，气候变化明显。冬春两季受西伯利亚和蒙古冷高压控制，气候干燥少雨，主导风向为西北风，风力大，风蚀外力地质作用极强。研究区北部的浑善达克沙地和科尔沁沙地生态环境极其脆弱，是北京西北部的主要沙尘暴源区。区域夏秋两季受太平洋副热带高压控制，东南风，风力弱，水汽供应少，气候炎热少雨。区域年平均气温在12.6℃以下，年降雨量200 ~ 750毫米，但降雨集中，降雨强度大，附加区地形坡度大，土质疏松，水蚀型外应力地质和重力侵蚀强烈，水土流失严重，易发生滑坡和泥石流。

2数据和预处理

2.1遥感数据

本文使用的遥感数据是美国EROS数据中心提供的MODIS影像。NDVI数据是2001到2004年6月6日合成的时间序列数据，有23个相位，空间分辨率为250m..地表温度(LST)是2002年8天的合成时间序列数据，有***46个相位，空间分辨率为1 km。

在MODIS数据处理中，利用MRT几何校正和镶嵌软件完成图像的几何校正和镶嵌。然后用最大合成法(MVC)对同一区域的植被指数、陆面温度等多时相数据进行合成预处理，即在J天内用最大像素值代替图像中的每个像素。这种处理的目的是减少大气云层、粒子、阴影、视角和太阳高度角的影响(布伦特，新罕布什尔，1986)。虽然最大合成过程(MVC)减少了大气中云和粒子的影响，但是云污染依然存在。然后，利用改进的最佳指数斜率提取方法对NDVI多时相相位差云进行处理。虽然MODIS的LST数据都是8天合成数据，但Ts数据的质量很差。为了解决数据不完全的问题，我们使用线性回归来模拟这些数据。地表温度在空间上有很强的相关性，同一区域的Ts在相邻时在空间上有某种相同的相关性，所以这种关系用线性关系拟合。使用相同大小的模板同时在恢复图像和参考图像上滑动。如果恢复图像中模板的中心值为零或异常，用最小二乘法求出两个模板中有效数据之间的线性回归系数，然后用该系数和参考图像模板的中心值求出新值，代替原来的零或异常值。

2.2其他辅助数据

辅助数据主要包括2002年北京市西北部土地利用现状图、北京市西北部1∶50000 DEM数据、北京市西北部降水和气温数据。根据北京西北地区气象站的资料，首先计算出各站的年平均积温和年平均降水量，然后用克里金插值法得到北京西北地区网格的年平均气温和年平均降水量分布图。

3研究方法

3.1土地覆盖分类

(1)选择适合MODIS数据分类的土地覆盖分类系统，本研究采用基于遥感数据的土地利用/土地覆盖分类系统[5]。这种分类体系最重要的特点是，对于不同的空间尺度和相应的遥感数据源，都有相应的分类，分类类型逐步细化。对于一级分类和二级分类，强调土地覆盖分类，即对于中低空间分辨率的遥感数据，以土地覆盖分类为主。

(2)利用主成分分析法对NDVI时间序列数据进行信息增强和压缩，消除各种干扰因素，提高分类精度。利用PCA变换，可以将原来12个月中的大部分有用NDVI信息压缩成少数几个第一主成分，并消除一些数据质量等原因造成的噪声。因此，PCA变换可以有效地保证分类精度不受损失。对实际结果的研究也表明，PCA在抑制噪声影响和保证分类精度方面发挥了重要作用[6]。

(3)结合LST数据、DEM数据和降雨温度数据，利用模糊K均值非监督分类方法[7]对研究区进行土地覆盖分类，分类处理后对分类误差明显的图斑进行校正，得到京西北地区土地覆盖分类图。

3.2多年来土地利用和植被覆盖变化分析

变化向量(CVA)分析是一种非常有潜力的植被对比分析方法，根据变化向量的强度和方向可以确定变化的区域和类型[8]。变化向量分析技术是将指标参数年度时间序列中的每个数据值作为时间序列空间中的一个点，将若干年的时间序列空间中的点连接起来，形成一个变化向量。变化向量的方向决定了变化的进度，向量的大小代表了变化的强度。

例如，如果用多年12个月的数据来分析变化向量，变化向量空间由每年12个变化监测指标的图像组成，那么全年的指标对应一个12维的时间向量:

土地信息技术创新与土地科技发展:2006年中国土地科学学会学术年会论文集。

P (i，X)表示像素I对应X年的向量，x (t)是像素I从时间t1到tn的指标值，n表示时间维度。向量的模‖P‖代表全年指标因子的累积，向量的方向是全年指标因子时间曲线形态的综合反映。

指标因子在任何两年的任何变化都会在这个12维空间中表现出来，可以用变化向量描述如下:

土地信息技术创新与土地科技发展:2006年中国土地科学学会学术年会论文集。

δP(I)是像素I从X年到Y年的变化向量。δP(I)包含像素I在(Y-X)年内每个时间维度的变化信息。变化向量的模‖δP(I)‖由Enclidean距离决定，Enclidean距离表示指标的变化强度。

土地信息技术创新与土地科技发展:2006年中国土地科学学会学术年会论文集。

当‖δp(I)超过一定阈值时，往往对应着植被覆盖类型从一种类型向另一种类型的变化。δP(I)的方向由一系列角度定义，决定了指标的变化过程。

根据影像和地面数据的直方图特征，可采用阈值分割法划分不同的矢量变化强度等级。

通过指示因子的累积值的变化率来判断向量变化类型。变化率定义如下:

土地信息技术创新与土地科技发展:2006年中国土地科学学会学术年会论文集。

4结论和讨论

4.1土地资源现状调查

通过主成分分析获得研究区2002年12个月的NDVI时间序列数据的前4个主成分，并基于MODIS 8天合成LST数据和1∶50000分辨率为1 km的DEM数据获得研究区的年平均LST数据。然后，结合降雨温度数据，采用模糊K-means无监督分类方法，获得北京西北部土地覆盖分类结果。然后对分类结果进行分类后处理，对分类误差明显的斑块进行修正，最终得到2002年北京市西北部土地覆盖分类图(图1)。

图1 2002年北京西北部土地覆盖分类图

从图1和图2可以看出，草地在北京西北部土地覆盖类型中所占比例最大，约占总面积的53%，内蒙古高原浑善达克沙地和丘陵区以及研究区东部科尔沁草原南缘集中连片；中西部的坝上高原低洼地区，河、湖、滩周边地区，阴山东部丘陵地区及周边地区也相对集中。草原总面积的60%以上分布在沙质甚至沙质干旱草原区。农业用地占研究区总面积的265，438+0%，主要分布在研究区西南部的高原和盆地，多呈带状沿河谷和河流冲积平原分布。林地占研究区总土地面积的65，438+03%，主要分布在大兴安岭、燕山、恒山和阴山地区。林地集中在研究区东部和西南部山区，大部分分布在山体上部。裸地占研究区总面积的8%，主要分布在北部的浑善达克沙地和东部的科尔沁沙地。在研究区域内，湿地、水域和建设用地所占比例最小。

图2 2002年北京西北部土地覆盖类型比例

4.2多年来土地利用和植被覆盖的变化

利用变化向量分析法监测了2001 ~ 2004年北京西北部NDVI的变化。所用数据为北京西北部2001至2004年的月最大NDVI时间序列值。首先计算NDVI的变化向量模，然后利用变化向量模的图像分割技术生成NDVI的变化强度。图像分割满足以下要求:①相似性原则，即同一区域的像素点应该相似；(2)不连续性原理，即从一个区域搜索到另一个区域时，像素的某些可变特性(如梯度)必须突然发生变化，从而确定边界。

4.2.1变化强度

NDVI的变化强度反映了植被覆盖的变化。综合考虑植被覆盖度变化向量模块的直方图、均值和方差，确定各个分割点，对变化向量模块进行分割，得到变化强度。

表1不同级别病媒变化强度的阈值

从图3和图4可以看出，2001到2004年的四年间，北京西北部大部分地区的土地利用/植被覆盖变化不大，生态系统基本保持平衡。无变化和低变化的面积占北京西北部总面积的92.3%。

未变化的面积最大，占总面积的53.7%，主要分布在北京西北部的赤峰市、敖汉旗、翁牛特旗、巴林右旗和阿鲁科尔沁旗，北部的锡林郭勒盟和四望子旗，表明近5年来三峡库区植被覆盖水平方向变化不大。

低变区占总面积的38.6%，主要分布在北京西北部的察哈尔右翼中旗、察哈尔右翼前旗、克什克腾旗、正镶白旗、正蓝旗、太仆寺旗交界处。

中变区面积相对较小，主要集中在北京西北部的凉城地区。

剧烈变化区主要集中在北京西北部西部、浑善达克沙地周围、克什克腾旗北部和四望子旗。

图3 2001-2004年北京西北部植被指数(NDVI)变化强度。

图4植被覆盖度矢量变化的强度-面积比

变更类型

上面计算了四年来北京西北部NDVI矢量变化的强度。矢量变化强度反映了2001到2004年北京西北部NDVI的变化程度，但无法判断这四年植被覆盖度是增加了还是减少了。因此，可以通过NDVI的变化强度和NDVI累积值的变化率来判断ND-VI向量的变化类型[9]。

阈值m是变化强度无变化和低变化的阈值。变化强度小于m的像素被认为是静止的。当变化强度大于m时，根据累积变化率确定变化类型。具体参数如下:

土地信息技术创新与土地科技发展:2006年中国土地科学学会学术年会论文集。

根据计算出的NDVI累积变化率，再考虑公式(5)，可以得到NDVI矢量变化类型图。

从图5和图6可以看出，从2001到2004年，北京西北部的植被覆盖度呈现出稳定增加的趋势。概括起来，变化特点如下:

图5 2001至2004年北京西北部植被覆盖指数(NDVI)。

图6植被覆盖度矢量变化的强度-面积比

(1)植被覆盖变化类型以固定型为主，占总面积的56.9%，主要分布在北京的东北部和西北部。

(2)增加型占比较大，占京西北地区的30.2%，主要分布在京西北地区的中南部。

(3)减少植被覆盖度有一定程度的减少，所占比例较小，主要分布在北京西北部以北的零星地区。

(4)波型占北京西北部面积的10.9%，主要分布在北京西北部的东北部。植被覆盖度的波动是一种正常的自然现象，是植被正常生长、长期气候变化和各种人类经济活动等自然作用的结果。

参考

[1]Tucker，C.J .，Townshend，J.R.G .和Goff，T.E .利用气象卫星数据进行大陆土地覆盖分类[J].理科，1984 (227):369~375

[2]Cihlar J .卫星大面积土地覆盖制图:现状与研究重点[J].国际遥感杂志，2000，21(6 & 7):1093 ~ 11113

[3]Townshend J.R.G，等.全球土地覆盖遥感分类:目前的能力和未来的可能性[J].环境遥感，1991 (35):243~255

[4]Lloyd D .利用短波植被指数图像进行陆地植被覆盖的物候分类[J].国际遥感杂志，1990，11(12):2269 ~ 2279

[5]J.Wang，T.He .发展基于遥感数据的中国土地利用/覆盖分类系统。SPIE会议录-遥感在环境监测、地理信息系统应用和地质学中的应用IV，2004，(5574):52~60

刘，，王静。基于PCA变换和神经网络分类的中国森林制图研究。长江流域资源与环境，2006，15 (1): 19 ~ 24。

刘，王静，等。基于MODIS数据的京西北地区土地覆盖分类研究。地理科学进展，2006，25 (2): 96 ~ 102

陈、、陈进、史。1983 ~ 1992中国陆地植被NDVI演变特征的变异向量分析。遥感学报，2002，6 (1): 12 ~ 18。

[9]中国测绘科学研究院。三峡库区生态环境监测技术研究。项目验收总结报告. 2005，10