基于MODIS地表温度和归一化植被指数的生态环境变化分析

贾虎军，杨武年，周丹，佘金星，许娟

(1.成都理工大学地学空间信息技术国土资源部重点实验室，成都 610059；2.成都理工大学遥感与GIS研究所，成都 610059)

1 引 言

植被是陆地生态系统的重要组成部分，其变化是研究全球环境变化的重要影响因子[1]。植被指数(Vegetation Index，VI)是基于地表反射率波段线性组合或原始波段比值，提取植被结构和功能信息的一种算法[2]，是指示大尺度植被覆盖的良好指标[3]。其中，Rouse等1974年提出的归一化植被指数(Normalized Difference Vegetation Index，NDVI)，可直接表征某一地区植被覆盖度，以及数据提取与处理过程中较强的抗干扰能力，已被广泛应用于植被变化分析[4-5]。中分辨率成像光谱仪(Moderate Resolution on Imaging Spectroradiometer，MODIS)卫星NDVI数据产品在地学综合研究中有了广泛的应用。MODIS NDVI在大范围植被时空变化[6-8]、森林火灾快速评估[9]、植被覆盖度动态遥感监测[10]等领域已做了深入的研究。NDVI和地表温度(Land Surface Temperature，LST)已被证明存在着一定的联系[11]，构成的特征空间蕴含着丰富、清晰的地表信息，可进行土地覆盖的变化和农业干旱评估等的研究[12-13]，是生态环境变化的重要指示因子。LST与NDVI之间的相关性在大洲和城市尺度已做了研究[14-16]，因此以中尺度与长年时间序列联合研究LST和NDVI具有重要意义。

本文采用普通线性回归分析方法研究2001年～2010年四川省MODIS LST和NDVI的变化趋势，揭示两者之间的内在联系。本文主要分为4个部分：①统计LST和NDVI 16天合成最大值的年平均值，分析两者的年际变化。②通过年平均值计算10年间LST和NDVI的变化率，以四川省21个市(州)分区统计回归斜率，比较每个市(州)正向和负向变化栅格数百分比。③以LST和NDVI回归斜率为数据源绘出二维散点图，研究两者的空间关系。④讨论在高程和坡度分级情况下，LST和NDVI的变化情况。

2 数据和方法

2.1 研究区概况

四川省位于中国西南部，地处长江上游，面积约48万km2，地理位置介于97°21′E～108°31′E，26°03′N～34°19′N之间。全省管辖18个地级市，3个自治州。全省地形西高东低，可分为川西高原和四川盆地两大部分[17]。因复杂地形和不同季风环流的交替影响，四川气候和植被复杂多样。东部盆地属亚热带季风气候，主要植被类型有常绿阔叶林和落叶阔叶林；西部高原在地形作用下，以垂直气候带为主，从南部山地到北部高原，由亚热带演变到亚寒带，垂直方向上由亚热带到永冻带的各种气候类型和植被类型[18]。

2.2 数据源

本研究使用的2001年～2010年MODIS NDVI数据是采用美国NASA(http：//lpdaac.usgs.gov)提供的MOD13Q1数据产品，该产品采用最大值合成法(constrained view angle maximum NDVI value compositing technique)，生成的250m、16天合成数据在一定程度上减少云、大气、太阳高度角和传感器等带来的影响，提高了数据的精度[7]。MODIS LST产品提供1km空间分辨率全球范围白天和夜间两次反演数据，该产品已被验证能够较好的拟合地表温度[19]。本文使用2001年至2010年LST产品的白天温度数据。在此次研究中共有MODIS LST和NDVI数据产品240期。

2.3 数据处理

利用Savitzky-Golay(S-G)滤波平滑工具对MODIS LST和NDVI时间序列数据进行平滑处理，可以有效去除LST和NDVI数据中的云、气溶胶及异常值的影响[2，20-22]。Savitzky-Golay滤波法是一种权重滑动平均滤波法，其权重取决于在一个滤波窗口范围内做多项式最小二乘拟合的多项式次数[20]。S-G滤波过程可由下式表示[21]：

(1)

其中，Y是原始序列NDVI值，Y*是重构序列NDVI值，Ci是第i个NDVI值滤波时的系数，N是指卷积的数目，m是窗口的宽度，j是指原NDVI数组的系数。

2.4 研究方法

首先由ArcGIS栅格计算器求得每年经Savitzky-Golay滤波平滑的LST和NDVI 16天合成最大值的年平均值，然后计算回归斜率。2001年～2010年间各栅格单元LST和NDVI变化趋势，根据线性回归斜率Xslope的计算公式来判读变化情况[23-24]。计算公式为：

(2)

其中，n为研究时间间隔，n=10；fi是第i年地表温度或NDVI均值；Xslope为线性回归斜率，若Xslope>0，则说明是增长趋势，即正向变化，若Xslope<0，则是下降趋势，即负向变化。

3 结果和分析

3.1 S-G滤波平滑

采用S-G平滑时，以窗口大小为9个像元数组与二次多项式重构LST和NDVI数据集，从而更好的描述LST和NDVI时序数据的总体变化趋势，研究结果更具有应用性。

3.2 LST和NDVI年际变化

根据图1得出，LST年平均值在14℃～17℃之间变化，呈小幅降低的趋势，并具有S型变化轨迹：2001年～2005年是LST下降时期；2006年～2010年是小幅度上升时期。NDVI 16天合成最大值的年平均值在0.58～0.63之间变化，呈现出与LST相反的变化趋势，且总体在小幅上升。从LST和NDVI年际变化来看，两者具有负相关性。

图1 2001年～2010年LST和NDVI 16天合成最大值年均值直方图

3.3 线性回归斜率分析

根据公式(1)计算四川省10年LST和NDVI线性回归斜率Xslope，由21个市(州)行政区划边界，分区统计各市(州)正向变化和负向变化栅格单元数，以及占栅格总数的百分比。表1中，正(n)或负(n)是正向变化或负向变化栅格数，正(%)或负(%)是正向变化或负向变化栅格数百分比。通过对表1的分析得出：各市(州)NDVI正向变化栅格数均大于负向变化，这与前文中分析得出的NDVI增长趋势是相符合的；而LST在攀枝花、绵阳市和广元市等8个市(州)正向变化小于负向变化。四川省自1999年开始实施退耕还林政策，到2010年全省森林覆盖率上升4%，NDVI值呈现正向变化。部分市(州)植被覆盖的恢复，对地表具有降温作用，LST出现负向变化。

根据表1绘出图2正、负向变化栅格数百分比差值曲线图，大部分市(州)LST和NDVI变化存在相反的趋势。即，当LST负向变化栅格百分比大于正向百分比时，NDVI负向变化百分比小于正向变化，反之亦然。其中，阿坝州、甘孜州和凉山州百分比差值较大的原因是其行政区划面积比较大。自贡市、泸州市、乐山市和宜宾市正向与负向变化栅格数百分比相等。

表1 线性回归斜率(Xslope)分区统计及百分比

图2 正向和负向变化栅格数百分比差值曲线

3.4 LST和NDVI的空间相关性

从图3散点图，可以看出各市(州)LST和NDVI大部分栅格单元回归斜率相关性较高，LST的变化率相对来说要大于NDVI。LST回归斜率在-0.2～0.2之间，而NDVI趋于-0.01～0.01之间。其中成都市和德阳较为分散，这是因为在城市快速发展中，大量原来植被覆盖区域变为城市用地，大量栅格单元NDVI变化率较大，接近于-0.02或0.02。城市的发展使不透水层面积扩大，而城市不透水层吸收太阳辐射加热地表和近地面大气，LST升高，变化率接近于0.4。大部分市(州)NDVI回归斜率Xslope>0，即植被覆盖度是增加的趋势。其中以广元市、广安市、达州市、巴中市和资阳市NDVI正向变化最为显著，LST则是负向变化较显著。值得注意的是泸州市和阿坝州的部分LST栅格单元斜率靠近0.4。

图3 LST和NDVI线性回归斜率(Xslope)散点图

3.5 LST和NDVI与高程和坡度的关系

LST和NDVI在很大程度上受到高程和地形因素的影响。本研究中通过自然间断点分割方法将高程分为：高程≤1000m、1000m<高程≤2000m、2000m<高程≤3000m、3000m<高程≤4000m和4000m≤高程五级；将坡度分为：坡度≤8°、8°<坡度≤16°、16°<坡度≤24°、24°<坡度≤32°和32°≤坡度五级。通过对图4的分析得出高程在4000m以下对于NDVI的影响较小，当高程大于4000m时NDVI值急剧减少。不同于NDVI的是，LST随着高程比较平缓的降低。不同年份在同一高程分级，LST和NDVI变化都比较小。由图5所示，当坡度小于等于8°时，NDVI值要大于坡度在9°～24°的值。坡度>24°时，因所占面积相对增加，NDVI值也在升高。LST随着坡度的增加而递减，且随着坡度增大降低速度变快。

图4 高程分级LST和NDVI值

图5 坡度分级LST和NDVI值

4 结束语

本文对利用Savitzky-Golay滤波技术平滑处理的四川省LST和NDVI值进行相关性分析，从年均值和线性回归斜率两个方面分析得出两者之间具有较高的负相关性。LST和NDVI年际变化具有相反的变化轨迹：若LST升高，NDVI降低；反之亦然。从各市(州)LST和NDVI二维散点图可以看出二者之间同样具有较大的相关性，即从空间角度分析LST和NDVI同样存在负相关性。综合来讲，以年为时间单位的中尺度LST和NDVI存在负相关性，这与前人以大洲和城市为尺度的研究结果相同。全省LST是轻微降低的趋势，而NDVI具有小幅的增加。四川省在全国率先实施退耕还林政策，森林覆盖提高，全省平均NDVI值有小幅增加，LST出现轻微降低的趋势。根据分级高程和坡度对LST和NDVI统计分析可知：4000m左右是高程对NDVI的影响结点，而LST随高程的增加而降低；LST和NDVI值在整体上随坡度增加而减少。

参考文献：

[1] 潘蔷，范文义，于海群，等.北京市植被指数变化与影响因素分析[J].北京林业大学学报，2012，34(2)：26-33.

[2] 李惠敏，刘洪斌，武伟.近10年重庆市归一化植被指数变化分析[J].地理科学，2010，30(1)：119-123.

[3] 李欣，王连喜，李琪，等.宁夏近25年植被指数变化及其与气候的关系[J].干旱区资源与环境，2011，25(9)：161-166.

[4] 田野，徐恒力，刘双，等.基于归一化植被指数变化分级的千山植被变化[J].中国水土保持，2010(01)：25-28.

[5] 樊风雷，邱茂慧.基于归一化植被指数的广州植被时间变化分析[J].广东农业科学，2012(06)：1-3.

[6] 林辉，熊育久，万玲凤，等.湖南省MODIS遥感植被指数的时空变化[J].应用生态学报，2007，18(3)：581-585.

[7] 王化儒，王天明，葛剑平.基于MODIS NDVI的东北虎分布区植被变化趋势[J].应用生态学报，2012，23(10)：2821-2828.

[8] 游浩妍，骆成凤，刘正军，等.基于MODIS植被指数估算青海湖流域植被覆盖度研究[J].遥感信息，2012，27(5)：55-60.

[9] 侍昊，王笑，薛建辉，等.基于MODIS NDVI时间序列的森林灾害快速评估方法[J].生态学报，2012，32(11)：3359-3367.

[10] 苗正红，刘志明，王宗明，等.基于MODIS NDVI的吉林省植被覆盖度动态遥感监测[J].遥感技术与应用，2010，25(3)：387-393.

[11] KAUFMANN R K，ZHOU L，MYNENI R B，et al.The effect of vegetation on surface temperature：A statistical analysis of NDVI and climate data[J].Geophysical Research Letters，2003，30(22)：43-55.

[12] 郭晓晓，冯海霞，冯仲科，等.基于MODIS温度和植被指数产品的山东省土地覆盖变化研究[J].遥感信息，2010(4)：58-62.

[13] SON N T，CHEN C F，CHEN C R，et al.Monitoring agricultural drought in the Lower Mekong Basin using MODIS NDVI and land surface temperature data[J].International Journal of Applied Earth Observation and Geoinformation，2012(18)：417-427.

[14] 梁保平，李艺，陈可宙.桂林市NDVI、地表温度的地物特征及相关性研究[J].遥感技术与应用，2012，27(3)：429-435.

[15] 马伟，赵珍梅，刘翔，等.植被指数与地表温度定量关系遥感分析[J].国土资源遥感2010(4)：108-112.

[16] JULIEN Y，SOBRINO J A，VERHOEF W.Changes in land surface temperature and DNVI values over Europe between 1982 and 1999[J].Remote Sensing of Environment，2006(103)：43-55.

[17] 郑泽忠，范东明，杨武年，等.四川省及成都地区土地利用/覆被变化建库与对比研究[J].测绘科学，2010，35(1)：143-145.

[18] 陈超，庞艳梅，潘学标，等.气候变化背景下四川省气候资源变化趋势分析[J].资源科学，2011，33(7)：1310-1316.

[19] 王之夏，南卓铜，赵林，等.MODIS地表温度产品在青藏高原冻土模拟中的适用性评价[J].冰川冻土，2011，33(1)：132-143.

[20] 吴文斌，杨鹏，唐华俊，等.两种NDVI时间序列数据拟合方法比较[J].农业工程学报，2009，25(11)：83-188.

[21] CHEN J，JÖNSSON P，TAMURA M，et al.A simple method for reconstructing a high-quality NDVI time-series data set based on the Savitzky-Golay filter[J].Remote Sensing of Environment，2004(91)：332-344.

[22] 边金虎，李爱农，宋孟强，等.MODIS植被指数时间序列Savitzky-Golay滤波算法重构[J].遥感学报，2010，14(4)：725-741.

[23] 贺映娜，白红英，高翔，等.基于NDVI的米仓山植被覆盖变化趋势分析[J].西北植物学报，2011，31(8)：1677-1682.

[24] 吴昌广，周志翔，肖文发，等.基于MODIS NDVI的三峡库区植被覆盖度动态监测[J].林业科学，2012，48(1)：22-28.