基于MODIS数据的羌塘国家自然保护区荒漠化遥感反演与评价

侯富强，李国明，阚瑷珂

（1.成都理工大学 信息科学与技术学院，成都 610059；2.成都理工大学 地球物理学院，成都 610059）

基于MODIS数据的羌塘国家自然保护区荒漠化遥感反演与评价

侯富强1，李国明2，阚瑷珂2

（1.成都理工大学 信息科学与技术学院，成都 610059；2.成都理工大学 地球物理学院，成都 610059）

基于MODIS数据采用因子分析法，选取MSAVI、NDVI、LST、TVDI四个指标反演羌塘国家自然保护区荒漠化动态变化过程，并进行分级评价。结果表明：①近十年间保护区荒漠化面积总体呈增加趋势，中度和极重度荒漠化面积呈减小趋势，轻度、重度荒漠化面积呈增加趋势；②若拉核心区、缓冲区东部无明显荒漠化现象，实验区与缓冲区西部荒漠化较为严重，美马错核心区、玛依核心区、缓冲区中部荒漠化程度较轻。

荒漠化；因子分析法；动态变化；分级评价

0 前言

荒漠化作为当今世界十大环境问题之首，倍受国际社会的关注［1］。我国是世界上荒漠化面积大、分布广、受其危害最严重的国家之一，西藏又是全国荒漠化较为严重的省份。藏北羌塘地区作为荒漠化监测和防治的典型区域，研究评价其荒漠化程度对高寒生态脆弱地区具有重要意义。作者在本文基于MODIS数据，采用因子分析法，选取了MSAVI、NDVI、LST、TVDI四个指标反演羌塘国家自然保护区荒漠化动态变化过程，并进行分级评价。

1 研究区概况

羌塘国家自然保护区（以下简称“保护区”，见下页图1）位于西藏西北部，昆仑山、可可西里山以南，冈底斯山以北，面积为29.8×104km2。范围涉及西藏那曲地区西部三县（安多、尼玛、班戈）和阿里地区北部三县（日土、改则、革吉）。根据《西藏自治区第三次荒漠化和沙化监测报告》资料可知以上区域是荒漠化较为严重的地区之一。

2 遥感反演指标选取

2.1 改进型土壤调整植被指数（MSAVI）

植被指数是根据植被反射波段的特性，计算反映地表植被生长状况、覆盖情况、生物量和植被种植特征的间接指标［2、3］。MSAVI作为改进型土壤调节植被指数，一般认为将土壤背景的影响减至最低，可增强对植被的敏感性。MSAVI指标采用植被指数数据产品MOD13Q1在2001年～2010年数据：MOD13Q1空间分辨率为250m×250m，时间分辨率为16d。保护区覆盖h24v05、h25v05共两个区域。

2.2 归一化植被指数（NDVI）

归一化植被指数（NDVI）能够敏感地反映出植被生物量、覆盖度，以及叶绿素含量等生物理化性质，在一定程度上反映了像元所对应区域的土地覆盖状况［4］，因而在土地覆盖变化研究中，常被用来描述土地覆盖的基本特征以及开展相关的定量研究。选用MOD13Q1的2001年～2010年数据，地面分辨率250m，每16天一次。

图1 羌塘国家自然保护区位置与区划Fig.1 The Qiangtang national nature reserve location and regionalization

2.3 陆地表面温度（LST）

对于植被茂密的地表，遥感反演所得到的地表温度，是指植被叶冠的表面温度。对于植被稀疏的地表，地表温度是地面、植被叶冠等温度的混合平均值［5］。采用MODIS数据的MOD11A2产品（2001年～2010年），时间分辨率8d，空间分辨率为1km的陆地表面温度产品，包含白天地表温度、夜间地表温度、31波段和32波段通道发射率等资料。

2.4 温度植被旱情指数（TVDI）

温度植被旱情指数为地表温度与NDVI的综合反映，被定义为基于NDVI－LST特征空间提取出的水份胁迫指标，较单独的地表温度指示度有所提高［6］。由于植被覆盖度与光谱植被指数存在一定关系，而植被覆盖度决定了传感器接收到土壤背景和植被冠层的可见光和热红外信息，从而影响了遥感影像获取的辐射温度［7］。土壤水份状况与表面温度之间不存在直接的关系，但土壤水份无疑是影响植被冠层温度的重要因素［8］。因此，一定植被覆盖条件下的冠层温度，能够间接反映土壤含水情况。

3 数据处理与分析

3.1 方法流程（见下页图2）

（1）针对MODIS数据的数据结构及通道特点，实现遥感数据的预处理、空间几何校正。

（2）针对MSAVI、NDVI、LST、TVDI四个指标进行荒漠化信息提取，对荒漠化过程进行反演。

（3）对各指标2010年与2001年荒漠化结果进行差值分析，完成基于各指标的荒漠化变化分布制图。

（4）分级评价羌塘自然保护区荒漠化空间分布和程度。

3.2 数据处理

3.2.1 MSAVI数据处理

使用ENVI提供的IDL接口，利用绿色植被在不同波段的光谱反射特征进行计算，由用户自定义函数计算得到。利用归纳法可推导得MSAV表达式，如公式（1）所示。

图2 方法流程Fig.2 Workflows of methods

式中 ρnir和ρred分别代表MODIS第1波段和第2波段的反射率。

反演结果见下页图3（a）。

3.2.2 NDVI数据处理

利用MRT（MODIS ReproJection Tool）对MOD13Q1－NDVI数据进行地理几何校正与重采样批处理，提取NDVI数据。同时，为了更加有效地消除云遮蔽、大气影响、观测中的几何关系等不利因素的影响，对保护区MOD13Q1－NDVI数据进行MVC（最大值合成法）处理，如公式（2），以获取每月NDVI最大值。

式中 NDVIi为第i月最大NDVI值；NDVIij为第i月第j个16－day合成NDVI值。然后，根据式（3），计算植被生长季累积NDVI。

式中 NDVIi为植被生长季累积NDVI；NDVI（i，k）为植被在第i年第k月最大NDVI值，i＝2001、2002、…、2010；k＝1、2、…、12。

反演结果下页图（b）。

3.2.3 LST数据处理

陆地表面温度是通过建立31通道、32通道亮温线性组合的分裂窗算法计算来获取。

MODIS数据共有八个热红外波段，根据分裂窗算法原理，MODIS数据选取第31、32两个波段来反演地表温度。采用覃志豪等［8］提出的，适用于MODIS数据的地表温度反演算法，该算法如公式（4）所示。

式中 LST是地表温度（K）；T31和T32分别是MODIS第31波段和32波段的亮度温度。根据这两个波段的图像DN值和普郎克公式来计算；A0A1、A2是分裂窗算法的参数，分别定义如下：

a31、b31、a32、b32是常量，可取：

a31＝－64.603 63；

b31＝0.440 817；

a32＝－68.725 75；

b32＝0.473 453。

其它中间参数分别计算如下：

其中 ti（q）是i（i＝31、32）波段视角为q的大气透过率。

反演结果如图3（c）所示（见下页）。

3.2.4 TVDI数据处理

温度植被干旱指数（TVDI）的定义为［6］式（10）。

式中 LST为陆地表面温度；Tmax＝a1＋b1＊NDVI，为某一NDVI对应的最高温度，即旱边；Tmin＝a2＋b2＊NDVI，为某一NDVI对应的最低温度，即湿边；a1、b1和a2、b2是旱边和湿边的线形拟合方程系数，由线性拟合得到；NDVI为植被指数。

反演结果如图3（d）所示（见下页）。

3.3 结果讨论

（1）MSAVI差值。值减小最为严重的地区在保护区的东南部，即实验区的东部、缓冲区的东北部，具体区域分布在班戈县，在保护区西南部表现为增加。

图3 2010年与2001年MSAVI/NDVI/LST/TVDI差值图Fig.3 The difference map of MSAVI/NDVI/LST/TVDI in 2001and 2010

（2）NDVI差值。负值越小，荒漠化程度越重正值越大，荒漠化程度越小。从下页图3中可以看出，其值减小最为严重的地区在保护区的东部，即实验区、缓冲区的东部、若拉核心区，主要区域分布在班戈县。

（3）LST差值。为便于研究，可将差值分为三类：①小于－1.00为温度下降区；②大于1.00为温度升高区；③之间的区域为微变区，其中升温集中的区域为尼玛县。

（4）TVDI差值。由于小于“0”的下降点非常少，所以不予考虑。在整个保护区范围内，TVDI指数都有不同程度上升，这表明旱情的增强。旱情明显增大的地区包括保护区的东北部、东南部和中西部地区，其它地区也有零星分布，其中数值较高的区域分布在班戈县境内。

4 分级评价

通过荒漠化程度样本图（国家林业局荒漠化检测中心绘制的2009年与2000年中国土地荒漠化图），作者与选取四个指标的MODIS数据进行GIS叠加分析，统计不同荒漠化程度的各指标值。通过对像元数值进行分析，确定保护区荒漠化分级评价体系（见下页表1）。

根据荒漠化分级评价体系计算各个指标的荒漠化分级面积，最后通过加权求平均值法求得2001年～2010年间保护区荒漠化程度分级统计面积及其变化趋势，如图4所示。

图4 羌塘自然保护区2001年～2010年荒漠化面积分级变化趋势Fig.4 The desertification area grading change trends o the reserve since 2001to 2010

保护区非荒漠化面积、中度荒漠化面积、极重度荒漠化面积呈减少趋势，轻度荒漠化面积、重度荒漠化面积呈增加趋势，且非荒漠化面积、极重度荒漠化面积减小较为明显，重度荒漠化面积增加较为明显。通过上述四种指标反演结果评价，保护区东部地区，即班戈县的大部份地区荒漠化程度较轻。该区域主要分布在湿润、亚湿润气候类型区降雨量较为丰富，但荒漠化程度却呈现较为严重趋势，同期该地区的年降水量也呈逐渐减少的趋势保护区西南部，即改则县中部、革吉县北部、日土县东部荒漠化情况较为严重，地处干旱、半干旱气候类型区，降雨较少，植被以草地为主，但荒漠化情况有所好转，年降水量增加。改则县北部，尼玛县的大部，班戈县的少部份变化较小，其气候比较湿润，降雨量均匀。影响荒漠化过程的除生态和气候等因素外，还来自于人类开发活动的干扰，过垦、过牧等，都是影响保护区荒漠化的可能因素。

表1 保护区荒漠化分级评价体系Tab.1 Desertification grading evaluation system of the reserve

5 结论

作者在本文中采用MODIS数据作为主要数据源，选取MSAVI、NDVI、LST、TVDI四项指标，对羌塘国家自然保护区近十年来荒漠化动态变化进行遥感反演，为高寒地区荒漠化遥感监测提供了可循的依据和方法。结果表明：

（1）非荒漠化面积在2001年～2005年间，以每年0.17×104km2的速度增加，在2005年～2010年间，以每年0.12×104km2的速度减少，总体呈减少趋势；中度和极重度荒漠化面积分别以每年0.03×104km2、0.05×104km2的速度减小；轻度、重度荒漠化面积分别以每年0.05×104km2、0.07 ×104km2的速度增加。

（2）若拉核心区、缓冲区的东部降雨量较为丰富，温度适中，没有呈现明显荒漠化现象；实验区的西部、缓冲区的西部降雨较少，温度较高，荒漠化较为严重；美马错核心区、玛依核心区、缓冲区中部荒漠化程度较轻。

［1］ 王莉雯，卫亚星，牛铮，等.基于GIS的沙漠化土地面积遥感分析监测［J］.地球信息科学，2007，9（3）：132

［2］ 王涛.干旱区绿洲化、荒漠化研究的进展与趋势［J］中国沙漠，2009，29（1）：1.

［3］ 张红，苏香玲.克拉玛依市土地荒漠化评价及动态监测［J］.安徽农学通报，2010，16（13）：235.

［4］ GU J，LI X，HUANG C，et al.A simplified data assimilation method for reconstructing time－serie modis ndvi data［J］.Advances in Space Research 2009，44（7）：501.

［5］ SHENG J F，WILSON J P，LEE S.Comparison o land surface temperature（LST）modeled with a spatially－distributed solar radiation model（SRAD）and remote sensing data［J］.Environmental Modelling ＆Software，2009，42（3）：436.

［6］ 李正国，王仰麟，周清波，等.基于TVDI方法的黄土高原地表水分时空分异研究［J］.中国农业资源与区划，2008，29（2）：13.

［7］ 周金星，张华，周泽福，等.大兴县荒漠化生态安全评价研究［J］.干旱地区农业研究，2009，27（1）：200.

［8］ 覃志豪，高懋芳，秦晓敏，等.农业旱灾监测中的地表温度遥感反演方法——以MODIS数据为例［J］.自然灾害学报，2005，14（4）：64.

book=45,ebook=45

1001—1749（2012）03—0353—05

X 144

A

10.3969/j.issn.1001－1749.2012.03.21

侯富强（1953－），男，工程师，从事计算机应用研究。

国家自然科学基金项目（40761018）

2012－04－17 改回日期：2012－04－19