植被指数与气象因子在像元尺度的多元线性回归研究

栾金凯, 刘登峰, 刘慧, 林木, 黄强

(1.西安理工大学 省部共建西北旱区生态水利国家重点实验室,陕西 西安 710048; 2.中国科学院地理科学与资源研究所,陆地水循环及地表过程重点实验室,北京 100101; 3.中国科学院大学,北京 100049; 4.中国水利水电科学研究院,北京 100038;5.中央财经大学 统计与数学学院,北京 100081)

植被在地球生态系统平衡、气候变化和水循环中起着协调的作用[1]。在流域的水文过程中,植被对蒸发、产流具有重要的影响。植被可在较短时间内反映气候变化和人类活动的影响[2]。作为植被生长及覆盖的最佳指示因子,归一化植被指数(Normalized Difference Vegetation Index,NDVI)是遥感数据的地表特征参量,反映植被分布特征及变化,是植被研究的基础数据之一[3-6]。

在有关植被指数影响因素的研究中,研究较多的是气温和降水,它们与植被指数的变化存在着密切的关系。ZHAO Xia等[7]对新疆1982—2003年的NDVI分析发现,新疆NDVI的增加与降水量和潜在蒸散发量的增加有关。FU Baihua等[8]分析了河岸植被NDVI与降水、气温、地表径流、洪水及地下水位之间的关系,发现最高气温是影响NDVI值最主要的因素,其次是前28 d的降水量。赵舒怡等[9]基于像元尺度,应用简单相关分析方法分析了2001—2013年植被指数与Palmer干旱指数(Palmer Drought Severity Index,PDSI)的相关关系。BIRTWISTLE Amy N等[10]研究发现,可以利用植被指数数据为干旱和半干旱无雨量站观测的偏远地区获取降水量。GEORGANOS Stefanos等[11]运用地理加权回归方法,对NDVI和降水数据进行了分析,发现可以使用降水数据对NDVI进行较好的预测。WEN Zhaofei等[12]基于像元尺度,分析了三峡库区1982—2011年的NDVI与降水、气温、太阳辐射及人类活动的相关关系。何云玲等[13]利用趋势分析和线性相关分析方法对云南省植被的月变化趋势、年际变化趋势进行了详细分析,并利用相关分析方法分析了NDVI与主要气候水热因子的关系。栾金凯等[14]应用复直线回归分析方法对榆林市近17年来的植被指数进行分析发现,在榆林市一半以上的区域内,人类活动对植被生长起到了促进作用;同时,还对汉江流域植被指数的影响因素进行了研究[15-16]。

目前,对NDVI影响因素的研究主要包括3个方面:①对研究区域的植被指数取面平均值,并与降水、气温等影响要素进行简单相关分析或回归分析;②在像元尺度上进行简单相关分析;③在研究区域内选取气象站点周围的NDVI值或面平均值,与气象因子和其他影响因子进行简单相关分析和多元回归分析。对NDVI进行像元尺度上的多元(3个及以上影响因素)回归分析的相关研究尚未见到。进行NDVI像元尺度的分析,可以准确地分析研究区域内每个地点的影响因素对小范围植被的影响,其空间连续性和异质性能够得到更好的体现和分析;利用多元回归分析,可以同时考虑多个气象因素对归一化植被指数的影响,可以比较准确地基于像元尺度预测未来植被的覆盖状况。因此,本研究以陕北地区的榆林市为研究区域,筛选出影响该区域归一化植被指数的4个气象因子,在像元尺度上对其归一化植被指数进行多元线性回归分析,分析回归系数的空间分布特征,评价回归分析得到的归一化植被指数值,为生态保护政策实施效果的监测提供参考依据。

1 研究区域概况

榆林市位于陕西省最北部,东临黄河与山西相望，西连宁夏、甘肃，北邻内蒙古，南接同省的延安市。经、纬度范围分别为107°28′E～111°15′E,36°57′N～39°35′N。地域东西长385 km,南北宽263 km,总土地面积43 578 km2。地貌大体以长城为界，北部为风沙草滩区，占总面积的42%，南部为黄土丘陵沟壑区，占总面积的58%[17-18]。榆林市的位置和地形如图1所示。

图1 榆林市位置和地形图

榆林市位于我国东部季风气候和西北干旱大陆性气候的过渡地带,也是毛乌素沙漠南缘与陕北黄土高原的过渡地带,这也决定了榆林市生态环境的脆弱性。榆林市位于“三北防护林体系工程”区域内,当前我国北方旱区森林覆盖已经接近气候容量的极限,极端干旱事件正在导致我国北方旱区普遍出现森林死亡现象,所以考虑到未来极端气候发生频率增加,榆林市不宜进一步增加森林覆盖。在这样的背景下,对榆林市植被的影响因素及未来发展趋势进行研究变得至关重要。

2 数据来源与处理

本文采用来自MODIS/Terra网站提供的NDVI遥感数据,数据集全称为MODIS/Terra Vegetation Indices 16-Day L3 Global 250 m SIN Grid V005,简称MOD13Q1。首先对下载的影像应用MRT(MODIS Reprojection Tool)软件进行批量裁剪、投影及转格式操作。本研究只对8月份的NDVI影像进行分析,MOD13Q1数据8月份有2幅影像,将得到的影像进行最大化合成(Maximum Value Composites,MVC)处理,最后每年得到一幅影像,所以2000—2017年共有18幅影像。

本文所用的8个气象数据(日照时数、平均本站气压、平均最低气温、平均最高气温、平均气温、平均风速、平均水汽压、平均相对湿度)来自中国气象数据共享服务网(http://data.cma.cn/),选取了榆林市内及其周围共17个气象站点的2000—2017年共18 a的气象数据,17个站点(如图1所示)分别为榆林、神木、定边、靖边、吴旗、横山、绥德、延安、河曲、兴县、离石、盐池、东胜、鄂托克旗、环县、五寨、隰县。利用反距离权重插值的方法来获得每一个像元的气象数据,插值出的气象数据的像元大小与NDVI像元大小相同且空间位置对应。

3 研究方法

3.1 相关性分析法

分别采用Pearson相关分析法和学生t检验法对NDVI与气象数据进行相关分析和显著性检验,其计算公式如下[19]:

(1)

(2)

3.2 多元线性回归分析法

本研究应用多元线性回归分析方法分析每个像元的NDVI与对应像元的多个气象因素之间的关系,多元线性回归分析的表达式为[19]:

θ=A0+A1x1+A2x2+…+Anxn+ζ。

(3)

式中:θ为NDVI值;A0、A1、…、An为回归系数;x1、x2、…、xn为气象因素,本文中n=4;ξ为残差,即实际NDVI值与线性回归分析结果的差值。

4 结果分析

4.1 气象因子的选取

本文选取了日照时数、平均本站气压、平均最低气温、平均最高气温、平均气温、平均风速、平均水汽压和平均相对湿度共8个气象因子。由于气象因子存在空间差异性,选取了榆林市内及其周围共17个气象站点资料,利用反距离权重插值的方法来获得每一个像元的气象数据,插值出的气象数据的像元大小与NDVI像元相同。考虑到气象因子对NDVI的影响存在滞后性,将榆林市2000—2017年8月份的NDVI分别与6月、7月、8月的各气象因子及6—7月、7—8月、6—8月各气象因子的平均值进行相关分析,相关系数见表1。由表1可知,与榆林市8月份NDVI相关性最大的各气象因子分别为7月日照时数、6—7月平均气压、7月平均最低气温、7月平均最高气温、7月平均气温、6月平均风速、7月平均水汽压和7月相对湿度。

表1 榆林市8月份NDVI与各时段气象因子间的相关系数

对这8个与NDVI相关性最大的气象因子进行t分布检验,结果见表2。由表2发现，6—7月平均气压、7月平均最高气温、7月平均气温、7月平均相对湿度通过了95%的置信度检验。因此,本文选择6—7月平均气压、7月平均最高气温、7月平均气温和7月平均相对湿度作为NDVI的主要影响因子。

表2 部分气象因子的t分布检验值

4.2 多元线性回归分析

选择上述17个气象站点2000—2017年6—7月平均气压、7月平均最高气温、7月平均气温、7月平均相对湿度,运用反距离权重法进行空间插值。插值后的栅格数据的像元大小和MODIS NDVI的像元大小相同,每幅影像为1 328×1 559的矩阵,去除周围不在榆林市境内的区域,共得到804 015个值,即榆林市共有804 015个像元。将2000—2014年共15 a的数据作为模拟数据用于建立模型,2015—2017年共3 a的数据用于验证模型。经过回归分析可得各像元的回归方程，其表达式的一般形式为：

θ=A0+A1x1+A2x2+A3x3+A4x4。

(4)

式中:θ为NDVI的回归分析值;x1为7月平均相对湿度;x2为6—7月平均气压;x3为7月平均气温;x4为7月平均最高气温;A0为回归常数;A1、A2、A3、A4分别代表7月平均相对湿度、6—7月平均气压、7月平均气温和7月平均最高气温的回归系数。

各回归参数(包括回归常数和各回归系数)的空间分布及频率分布分别如图2和图3所示。由图2可知：回归常数在大部分区域为正值,在榆阳区和神木市部分区域为负值,西部的定边县和东南部的部分区域为0;7月平均相对湿度的回归系数在榆林大部分区域为正值,说明在大部分区域,制约榆林植被生长的主要因素是水。从空间分布来看,6—7月平均气压与7月平均相对湿度对植被的影响基本呈现相反的趋势,7月平均气温与7月平均相对湿度的回归系数也呈现出相反的趋势,而7月平均最高气温的回归系数与7月平均气温的回归系数呈现出相反的趋势;不同的影响因子对同一地区的响应不同,不同区域对同一影响因子的响应也不同。

图2 各回归参数的空间分布图

图3 各回归参数的频率分布图

表3为各回归参数(804 015个值)的基本统计特征值,包括平均值和标准差。

由表3可知,各影响因素的系数分布分散,说明不同地点上影响因素的作用有较大差异。

表3 各回归参数的基本特征统计值

4.3 NDVI模拟值的评价

榆林市2015—2017年各年NDVI模拟值与观测值的空间分布如图4所示。分析图4发现：榆林市各年NDVI的模拟值与观测值的空间分布特征基本一致,毛乌素沙漠边缘NDVI值较低,榆林东部和东南部NDVI值较高;整体的模拟效果较好,其中2016年的模拟效果最好,2017年的模拟效果较差。

图4 榆林市2015—2017年模拟NDVI与实际NDVI的空间分布图

为了定量评价逐个像元的模拟效果,将实际NDVI值减去模拟NDVI值的结果记为残差。残差可视为是人类活动对归一化植被指数的影响。如果残差为负值,说明人类活动对该区域的植被起到了抑制作用,可能出现了滥砍滥伐、植被破坏或者植被死亡等现象。如果残差为正值,说明人类活动对植被的生长起到了积极的促进作用。残差的绝对值与实际NDVI结果的比值记为相对误差,是评价模拟效果的关键指标。相对误差越大,说明人类活动对植被的影响越大;反之越小。

榆林市2015—2017年各年NDVI模拟值的残差和相对误差的空间分布如图5所示。从图5可以看出：2015年残差为负值的区域主要位于靖边县、横山区、子洲县、神木市的北部及府谷县的南部,说明在这些区域存在着植被破坏或死亡的现象;其余区域残差值均为正值,说明在这些区域内,人类活动对植被生长起到了积极的促进作用,这些地区实施退耕还林、荒山荒地造林、封山育林等措施的效果较好。2016年残差为负值的区域主要分布在神木市的北部和府谷县的南部。2017年残差为负值的区域主要位于靖边县的西部,佳县、米脂县和子洲县的部分区域。

图5 榆林市2015—2017年各年NDVI模拟值的残差和相对误差分布图

榆林市2015—2017年各年NDVI模拟值的残差及相对误差的统计情况分别见表4和表5。由表4和表5可知：对于榆林市NDVI的模拟值,2015年的残差在[-0.1,0.1)上的像元比例为66.50%,相对误差在25%以内的像元比例为68.68%;2016年的残差在[-0.1,0.1)上的像元比例为77.62%,相对误差在25%以内的像元比例达到了88.93%;2017年的残差在[-0.1,0.1)上的像元比例为72.82%,相对误差在25%以内的像元比例为89%。由此可以看出：2016年的模拟效果最好,该年人类活动对植被的影响最小;2015年的模拟结果最差,该年人类活动对植被的影响最为剧烈。

表4 榆林市2015—2017年NDVI模拟值的残差统计

表5 榆林市2015—2017年NDVI模拟值的相对误差统计

由表4和表5还可以得知,榆林市2015—2017年NDVI模拟值的残差为正值的像元比例分别为60.14%、79.31%和85.73%,说明最近几年人类活动对植被生长和更新起到促进作用的区域占到了很大部分。查阅榆林市林业局网站发现[20],榆林市2015年总造林面积45万亩,2016年和2017年分别造林73.3万亩和70万亩,这也导致2016和2017年残差为正值的区域面积大于2015年的。另外,2016年横山区造林面积15 300万亩,所以该区NDVI模拟值的残差由2015年的负值变为了2016年的正值。2016年春季,府谷县存在严重的偷牧乱牧现象,所以2016年该区NDVI模拟值的残差值小于0。2017年7月,绥德县和子洲县受到暴雨和洪灾的影响,当地林业受损严重,导致2017年这两个县NDVI模拟值的残差值呈现负值。由此可见,应用本文所采用的方法,可以作为监测封山育林、退耕还林、退牧还草等政策的落实和实施效果的参考;在现阶段人类活动变化幅度不大的情况下,还可以根据未来的气象条件预测未来植被理论上的生长状况。

对榆林市2015—2017年NDVI的实测值与模拟值进行逐像元相关系数分析,并作实际值与模拟值的散点图,横坐标为榆林市804 015个像元的NDVI实测值,纵坐标为其模拟值,如图6所示。由图6可知：2017年的模拟效果最好,相关系数达到了0.851 0;2015年的相关系数为0.616 6;2016年的相关系数为0.803 3。

图6 榆林市2015—2017年NDVI实测值与模拟值的散点图

5 结论

本文利用MODIS/Terra NDVI时间序列数据,基于像元尺度对陕西省榆林市2000—2017年生长季(8月份)的归一化植被指数(NDVI)应用多元线性回归方法进行影响因素分析研究,主要得到以下结论:

1)与榆林市8月份NDVI相关程度最大的因子为:6—7月平均气压、7月平均最高气温、7月平均气温和7月平均相对湿度,可视为榆林市8月份NDVI的主要影响因子。

2)不同的影响因子对同一地区的响应不同,不同区域对同一影响因子的响应也不同。7月平均相对湿度与6—7月平均气压对植被的影响基本呈现相反的趋势。7月平均气温与7月平均相对湿度的回归系数呈现出相反的趋势,7月平均最高气温的回归系数与7月平均气温的回归系数也呈现出相反的趋势。

3)对比NDVI的模拟值和实测值发现:整体的模拟效果较好,2015年NDVI模拟值的残差为负值的区域主要位于靖边县、横山区、子洲县、神木市的北部及府谷县的南部,2016年NDVI模拟值的残差为负值的区域主要分布在神木市的北部和府谷县的南部,2017年NDVI模拟值的残差值为负值的区域主要位于靖边县的西部,佳县、米脂县和子洲县的部分区域。残差为正值的区域,说明人类活动对该区域内植被的生长起到了促进作用,这些地区的封山育林、退耕还林、退牧还草等措施的实施效果较好。残差值为负值的区域,说明人类活动对该区域内植被的生长起到了抑制作用,可能存在乱砍乱伐、植被死亡等现象。应用这种方法,在现阶段人类活动变化幅度不大的情况下,可以根据未来的气象条件预测未来植被理论上的生长状况,可以为区域生态环境保护效果的评价提供参考。