贵州乌江流域NDVI降尺度适应性评价

2021-01-20杨江州周旭熊军周志凯韦小茶

遥感信息 2020年6期

杨江州，周旭，熊军，周志凯，韦小茶

(1.贵州省地质矿产勘查开发局 117地质大队，贵阳 550018；2.贵州师范大学地理与环境科学学院，贵阳 550025)

0 引言

植被是陆地生态系统的主体，具有广泛性、复杂性、丰富性的物种资源，是地表水循环、生物循环的重要途径[1-2]。随着遥感技术的发展，其观测面积广、信息量大、分辨率高、时间序列长等优势，在植被监测中得到广泛应用[3-5]。归一化植被指数(normalized difference vegetation index，NDVI)是目前公认的表征植被变化的最有效参数之一，已广泛用于植被生长及其动态变化研究中[5]。王茜等[6]应用NDVI数据研究我国30年NDVI时空变化，并对气候因素影响进行分析；蒙吉军等[7]基于NDVI、NPP数据研究西南喀斯特植被变化对气候的响应；朱林富等[8]对山地城市进行植被变化时空分异特征分析；张继等[9]基于NDVI对生态工程建设背景下贵州高原植被变化及因素分析；王冰等[10]基于NDVI数据分析贵州喀斯特地区植被覆盖变化趋势。在降尺度研究中，统计降尺度方法得到广泛运用。如Huth[11]运用神经网络与线性回归方法对欧洲日气温进行降尺度研究，发现较为简单的线性方法也可以取得较好的降尺度效果；Immerzeel等[12]通过建立TRMM降水数据与NDVI 关系，对伊比利亚半岛的TRMM年降水数据进行降尺度研究，将RMM 数据的空间分辨率从 0.25°提高至 1 km。之后Duan等[13]、李净等[14]、范雪薇等[15]在此基础上引入地形因子进行多元回归，实现了TRMM降水数据降尺度，且能够体现更细致的降水分布特征。然而在植被变化的研究中，高精度、小尺度的时间序列均较短，多为2000年至今的研究，对植被降尺度研究较为缺乏，且1982—2000年间NDVI产品多是8 km低分辨率数据。表明利用该尺度数据揭示植被变化仅能反映大范围的总体特征，小尺度探究植被变化存在不足。

贵州乌江流域地理环境空间异质性大，高分辨率遥感数据揭示中小流域NDVI变化及其驱动因素尤其重要[16-17]。乌江流域地处中国西南喀斯特腹地，岩溶发育典型，土壤瘠薄，水资源赋存条件极弱[18-19]。流域特殊的地质背景，再加上人类活动影响，生态环境不断退化，使其演化为典型生态脆弱区[20-21]。因此，研究喀斯特流域NDVI长时间变化趋势，有助于揭示生态环境演化特征，减弱不合理人类活动对其的影响，对山地生态文明建设具有积极意义。本文基于空间统计降尺度方法，探究喀斯特流域NDVI的降尺度与时空变化。研究成果对低分辨率NDVI数据降尺度获取,喀斯特地区石漠化治理与生态建设具有重要指导意义。

1 研究区概况、数据与方法

1.1 研究区概况

贵州乌江流域位于中国西南腹地，自西向东横贯整个贵州省北部，向西与云南省接壤，向北与四川省和重庆市接壤。流域干流全长847 km，是贵州省最大流域，支流较多，呈羽状水系分布。流域属亚热带季风性湿润气候，气候温和，降水丰沛，雨热同季，年平均气温15 ℃，年平均降水量约1 200 mm。流域地质构造上主要属于扬子准地台中的黔北台隆，在地域独特的水热条件的驱动下，碳酸盐岩被流水侵蚀，使流域内岩溶地貌发育强烈，75.6%的地区是石灰岩、白云岩等碳酸盐岩[22]。由于碳酸盐岩抗侵蚀能力较强，母质造土能力差，成土过程缓慢，不利于土壤的形成与植被生长，使其基岩裸露率高，石漠化敏感性强。地貌类型主要分为上游的岩溶峡谷地貌，中游的岩溶高原地貌，下游的岩溶槽谷地貌，平均海拔约1 100 m。

1.2 数据来源

GIMMS NDVI 数据来自NOAA网站(https：//www.noaa.gov)网站下载，空间分辨率为8 km，时间分辨率为15 d，时间序列为1982—2015年；MODIS NDVI数据来自美国NASA网站(https：//modis.gsfc.nasa.gov)网站下载MODIS13 Q1，空间分辨率为250 m，时间分辨率为16 d，时间序列为2000—2016年。

1.3 研究方法

本文采用的NDVI数据是基于GIMMS和MODIS13 Q1 2种遥感产品获取，都源于美国国家航空航天局(National Aeronautics and Space Administration，NASA)的一系列遥感卫星。GIMMS NDVI数据集是使用时间较长的植被数据，该数据空间分辨率为8 km，时间分辨率是15 d，下载1982—2015年数据，总共816幅影像，利用Matlab编程语言转换GIMMS NDVI3g.v1原始产品为Geottif格式，采用WGS_1984地理坐标和Albers Equal-Area Conic投影坐标。MODIS NDVI数据是植被数据中具有较高分辨率且应用广泛的数据，空间分辨率为250 m，时间分辨率为16 d，下载2000—2016年数据，总共408幅影像，利用MRT(MODIS reprojection tools)软件对数据进行投影转换、数据拼接等预处理，采用与GIMMS NDVI数据一致的地理坐标和投影坐标。再运用ENVI/ArcGIS对数据集进行裁剪、拼接、几何校正等预处理工作，同时剔除NDVI产品的无效填充值。采用最大值合成法(maximum value composites，MVC)将每月中2期NDVI影像提取最大值，作为每月的NDVI值，此方法能进一步消除大气、卫星轨道飘移、太阳高度角等的干扰，在地形复杂的山地具有较好的效果。其计算如式(1)所示。

NDVIi=Max(NDVIij)

(1)

式中：NDVI为归一化植被指数；NDVIi指第i个月或者第i年的NDVI；NDVIij指第i年的第j个半月的NDVI或第i年第j月的NDVI。

回归模式法是应用最早和最广泛的统计降尺度技术，是通过建立预报变量与预测因子之间线性或非线性关系的一种概念性方式[11]。应用IDL(interactive data langguage)编写线性回归方程代码，实现NDVI数据降尺度，获取更早时间序列及较高空间分辨率的数据。其原理利用2001—2015年的GIMMS NDVI与MODIS NDVI数据进行15年逐像元回归分析，且该GIMMS NDVI数据的空间分辨率都重采样为250 m分辨率，与MODIS NDVI分辨率一致。得出2001—2015年的GIMMS NDVI与MODIS NDVI数据像元回归方程的参数后，基于1982—2000年GIMMS NDVI数据，反推出各像元上的NDVI值，从而得到1982—2000年250 m的NDVI数据，具体流程如图1所示。式(2)为线性回归方程。

GNi=a+b×MNi+ci

(2)

式中：GNi代表第i月、年的GIMMS NDVI；MNi代表第i月、年的MODIS NDVI；a和b为常数；c为回归方程的残差。

图1 降尺度流程图

采用平均绝对误差(M1)、平均相对误差(M2)、均方根误差(M3)、Theil不等系数(M4)探究GIMMS NDVI数据空间降尺度得到的NDVI和MODIS NDVI数据的精度，对其做时间序列与空间的可靠性分析(式(3)～式(6))[23-24]。

(3)

(4)

(5)

(6)

式中：ai实际值；bi为模拟值；n为检验点的个数；M4介于0和1之间，M4越趋近0，则代表实际值与模拟值误差越小，模型模拟的精度越高。

利用线性趋势法进行NDVI时间趋势的分析，模拟出每个栅格的变化趋势，以单个像元时间变化特征综合反映出整个区域的时空演变格局和空间变化规律。随着时间序列的改变，NDVI表现为序列整体的上升或下降变化趋势，呈现出空间分布格局的变化和在某时刻出现突变或转折，NDVI变量可以看作是时间的一元线性回归，利用最小二乘法逐像元计算NDVI变化斜率如式(7)所示[25-26]。

(7)

式中：B为某像元的趋势线斜率；t为年份；n为研究年段35(1982—2016年)。当B值为正值时，表现为随时间t变化NDVI呈上升趋势；反之，当B为负值时，随时间t变化NDVI呈下降趋势。其绝对值越大表示研究区域NDVI上升或下降的趋势越显著。

2 统计降尺度模型结果验证

2.1 基于降尺度数据的可靠性分析

由于GIMMS与MODIS数据在2000年开始重叠，以2000年GIMMS NDVI降尺度为例，更能说明降尺度的效果特征，如图2所示。二者对比，整体来看，分辨率为8 km的GIMMS NDVI降尺度后得到分辨率250 m的NDVI影像更加清晰，更有利于揭示流域NDVI变化；虽然在GIMMS NDVI中，乌江流域西部和中南部的低值在统计降尺度后有一定的变化，但整体的趋势都呈现西北低、东南高的变化趋势。表明降尺度的NDVI整体分布状况与GIMMS NDVI在空间上有着较好的一致性。

图2 2000年贵州乌江流域GIMMS NDVI降尺度效果图

2.2 降尺度数据的时间验证

利用1982—2000年的GIMMS NDVI数据和降尺度NDVI数据的年际平均值统计尺度和多年各月平均值统计尺度进行时间相关性检验。年际平均值统计尺度相关性验证为1982—2000年全省平均的年时间序列，有19个时间点；多年各月平均值统计尺度相关性验证为1982—2000年多年各月平均值组成的站点序列，有12个时间点。

以年尺度时间序列来看，运用1982—2000年的GIMMS NDVI数据和降尺度NDVI数据做误差检验分析(表1)，结果表明：从实际误差(M)来看，误差值的正数值与负数值的个数各在一半，且多年总和为0.007 11，表明GIMMS数据降尺度得到NDVI数据与实际波动较为合理，误差较小。在绝对误差与相对误差中，1993年绝对误差(M1)与相对误差(M2)达到最高，分别为0.013 02、0.022 95；1992年绝对误差(M1)与相对误差(M2)值最低，分别为0.000 22、0.000 37。从数值趋势来看，90年代与80年代相比，具有上升的趋势。19年间，二者均方根误差(M3)为0.006 52，Theil不等系数(M4)为0.005 45，虽各年份的误差均有波动，但整体较小。

表1 年际平均值误差检验

从月尺度时间序列来看，运用1982—2000年各月均值的GIMMS NDVI数据和降尺度NDVI数据做误差检验分析(表2)，结果表明：实际误差(M)中，有7个正值，5个负值，大致各在一半，同样表明降尺度数据的合理性。从绝对误差(M1)与相对误差(M2)来看，6月的多年各月平均值绝对误差(M1)与相对误差(M2)达到最高，分别为0.084 32、0.126 32；原因可能在于6月是植被茂盛时期，不同植被每年长势不一样，而降尺度数据的NDVI被多年数值回归后导致误差偏大；9月绝对误差(M1)与相对误差(M2)值最低，分别为0.000 70、0.001 00。多年各月平均值均方根误差(M3)为0.029 56，Theil不等系数(M4)为0.029 45，相对于年尺度变化可发现，误差有所增大，但整体降尺度的月份趋势与GIMMS的月份趋势大致一样，可以反映月份序列的变化。

表2 多年(1982—2000年)各月平均值误差检验

图3 年际平均值统计尺度和多年各月平均值统计尺度拟合分析

根据GIMMS NDVI数据和降尺度NDVI数据做年际平均值统计尺度和多年各月平均值统计尺度相关性分析(图3)。年际平均值统计尺度相关分析的线性方程为y=0.635 6x+0.215，R2=0.998 9；多年各月平均值统计尺度相关分析的线性方程为y=1.238 6x-0.139 8，R2=0.949 6，都且在0.05显著性水平上显著相关。表明降尺度得到的NDVI数据与GIMMS NDVI数据拟合度较好。

综上所述，年际平均值统计尺度和多年各月平均值统计尺度都具有拟合精度较高，良好地反映时间变化趋势，说明降尺度得到的NDVI数据在时间序列上具有一定的代表性。但结合表1和表2比较来看，年际平均值统计尺度模拟更优。原因可能在于年际平均值统计尺度模型中，2000—2015年进行回归数据较少，一年1个数据，总共15个数据；而多年月均值统计尺度，一年模拟12个数据，总共180数据。表明多数据线性回归后，得到的参数，模拟NDVI具有均一性，而各年月气候环境的不确定性，使之降尺度得到的NDVI值误差精度增大，但整体误差值较小，误差精度满足NDVI时间变化的研究需要。

2.3 降尺度数据的空间验证

以2000年MODIS NDVI数据与2000年降尺度NDVI数据进行空间上的相关性检验，进行空间误差检验分析(图4)。从图4可知，2000年降尺度的NDVI值相对于MODIS数据来说，最大值从0.82变为0.84，相差0.02，整体值有所增加。但整体均值相差不大，差值为0.003 96，变化趋势保持较好的一致性。

图4 2000年贵州乌江流域MODIS NDVI与降尺度NDVI分布图

根据以上原理方法，利用ENVI对栅格数据进行波段运算，求出降尺度NDVI值与MODIS NDVI值的平均绝对误差(M1)、平均相对误差(M2)、均方根误差(M3)、Theil不等系数(M4)，进行统计分析得出表3，再运用ArcGIS进行数据处理得出空间图(图4)。由于NDVI绝对值处于0～1之间，采用理论值为小于等于5%的误差，视为精度效果好；理论值为5%～10%之间的误差，视为精度效果较好；大于等于10%的误差，说明误差较大。因此，本文按照实际情况误差值的大小分为3级：小于等于5%，精度良好；5%～10%，精度较好；大于等于10%，精度较差。

表3 2000年MODIS NDVI值与降尺度NDVI值误差占比 %

结合表3与图5来看，可知平均绝对误差(M1)、均方根误差(M3)、Theil不等系数(M4)空间占比小于等于5%误差值的大部分集中于中部与东北部，误差占比分别为80.82%、80.82%、86.20%；5%～10%误差值集中于西部，误差占比分别为18.11%、18.11%、13.02%；较少出现大于等于10%误差值，占比分别为1.07%、1.07%、0.80%。虽然平均相对误差(M2)大于等于10%的误差值有所增加，占比为15.50%，集中分布西部，在中部和东部都有零散的分布，但整体小于等于10%的占比为85.85%，整体趋势精度较好。表明降尺度的NDVI值可以良好地反映在空间的动态趋势。

图5 2000年贵州乌江流域MODIS NDVI值与降尺度NDVI值误差空间检验

3 植被时空变化特征及其趋势变化分析

3.1 植被年际时空变化特征

以1982—2000年降尺度NDVI数据与2001—2016年MODIS NDVI数据相结合，统计出乌江流域1982—2016年NDVI的年际变化，如图6所示。从整体趋势来看，多年的线性方程为y=0.001 4x-2.185 6，R2=0.321 6，贵州乌江河流域NDVI变化呈明显增加波动变化趋势。根据乌江流域NDVI的年际变化，可划分为2个阶段：第1个阶段为1982—2000年，流域NDVI值从0.590 1下降到0.564 2，线性方程为y=-9E-0.5x+0.764 3，R2=0.001 6，有向下的波动变化趋势，但并不明显；第2个阶段为2000—2016年，流域NDVI值从0.564 2上升到0.673 1，线性方程为y=0.006x-11.508，R2=0.827 7，NDVI值有明显的增加趋势。其原因在于，20世纪80年代处于改革开发初期，乌江流域居民对植被保护的意识不高，再加上人类活动对植被的乱砍滥伐，造成第一个阶段不明显的下降波动趋势。90年代末期，国家加大西南喀斯特地区的石漠化治理，使乌江流域植被逐步恢复改善，因此形成第2阶段明显增加的波动趋势。

图6 1982—2016年贵州乌江流域植被年际变化

应用ArcGIS对每年的NDVI影像进行栅格运算，最后得出多年平均NDVI空间分布，如图7所示。贵州乌江河流域逐年NDVI值具有明显的空间分布差异，由于受水平地带性、气候和地形地貌等因素的综合影响，呈现出东高西低的分布格局，NDVI值为0.06～0.80。NDVI低值区域主要分布于流域中西部，从西部区域看，形成低值区域原因为该区域处于贵州地势较高山区，热量较低，年降雨量偏少，不利于植被生长；从中部区域来看，低值区域主要集中于贵阳市与遵义市周边，其中贵阳市周边的低值区域最大，主要原因是城市的开发增强，土地利用格局发生深刻改变，影响了植被生长。这表明人类活动频繁区域对植被生长起到抑制作用，因此在开发建成区中，应加强城市建成区的绿化，改善人居生态环境。

图7 1982—2016年贵州乌江流域NDVI多年平均空间分布

3.2 植被趋势变化分析

基于利用线性趋势法，在像元尺度上分析了贵州乌江河流域NDVI年际变化趋势(图8)。按照实际情况趋势线斜率值的大小分为5级：B<-0.001为显著下降；-0.001≤B<-0.000 5为轻微下降；-0.000 5≤B≤0.000 5为基本不变；0.000 50.001为显著上升。从1982—1999年趋势来看，基本不变的趋势占主体，占总面积的43%，而上升趋势与下降趋势呈零散分布状，其中上升趋势占总面积的42%，下降趋势占总面积的15%，整体呈增加趋势。原因在于20世纪90年代石漠化治理在21世纪初期取得了初步成效，生态坏境得到了改善，植被覆盖度增加。

图8 贵州乌江流域NDVI年际变化趋势

4 讨论与结论

4.1 讨论

乌江流域地处贵州喀斯特高原，多云雨雾天气，覆盖研究区的Landsat等高分辨率遥感数据云量较多，数据质量较差，难以得到时间序列连续数据，而SPOT、AVHRR等数据满足长时间序列，但是空间分辨率极低，再加上喀斯特山区地形特征差异显著，可进入性极差，地表实测数据等基础资料相对缺乏。本文基于统计降尺度模型，将1982—2000年GIMMS NDVI数据的8 km分辨率提升到250 m，且与2000—2016年MODIS NDVI数据(分辨率为250 m)相结合，从而得到1982—2016年250 m的长时间序列高分辨率NDVI数据。通过该方法获得研究区长时间序列高分辨率的NDVI数据，高分辨率的NDVI数据更能提高喀斯特小空间尺度的植被变化参数特征，有利于揭示喀斯特山区流域植被变化及其驱动因素[14]。因此，基于NDVI数据建立2种各像元间的统计关系，实现数据统计降尺度，从而延长插补时间序列获取高分辨率长时间序列数据的方法。该方法也可对今后在难以获取长时间序列数据或小空间尺度需求更高分辨率的相关研究提供参考。

王秀春等[27]基于1998—2006年SPOT生长季(4—9月)逐旬最大合成NDVI数据，分析了乌江流域近10年的植被覆盖动态变化，结果数据显示 NDVI年际变化呈增加趋势，年增长率为0.009/a，流域内地级市及省会城市呈减少的趋势。在文献[9-10]基于NDVI数据对贵州喀斯特地区植被覆盖变化研究中，均发现贵州喀斯特地区总体处于上升趋势，降水对植被变化影响减弱。相对于前人的研究，本文以乌江流域为研究区，得出1982—2016年研究区的NDVI年际变化同样呈增加趋势，年增长率为0.006/a；流域植被年际趋势分布变化中，人类活动较大的地方植被变化呈减少的趋势；从整体来看，本研究的结果与前人研究基本一致。

在1982—2016年贵州乌江流域植被年际变化中(图6)，NDVI在2000年出现明显下降，原因可能为该年降水量偏小，年均温度较低，不利植被生长。该年现象与文献[28-29]研究NDVI变化趋势基本一致。但本文方法只是数据的交叉检验，没有辅助数据的佐证，缺乏一定的物理机制，在未来模型可以通过加入第3个辅助数据来进一步的探讨。同时也存在一定的不确定性，如 NDVI 在高植被覆盖区易出现饱和的问题；不同传递函数与降尺度结果精度的关系等问题。