基于地理加权回归模型的能源“金三角”地区植被时空演变及主导因素分析
2018-08-29李晶晶闫庆武胡苗苗
李晶晶, 闫庆武, 胡苗苗
(中国矿业大学环境与测绘学院, 江苏 徐州 221116)
植被可在全球气候变化研究中充当“指示器”[1],其覆盖变化能准确反映区域性生态环境状况,因此植被的退化和恢复问题日益受到国内外学术界的广泛关注[2]。归一化植被指数(normalized differential vegetation index,NDVI)是植被长势及其覆盖程度的最佳指示因子[3-5],并与地表植被覆盖率呈正比;对于同一种植被,NDVI值越大,表明地表植被覆盖率越大,植被长势越好。目前已有采用遥感技术研究不同类型区域植被动态变化的报道[6-8],已有学者针对干旱半干旱地区展开相关研究,如我国的三江源区域[9-10]和青藏高原地区[11]等。已有研究主要采用传统的统计回归及相关性分析方法探索植被变化及其与影响因子之间的相关性,进而研究这一相关性在不同区域、不同植被类型中的差异。但是,地理环境具有空间异质性,前人研究表明即使在相似的覆被条件或植被类型中,这种相关性在空间上也具有非平稳性[12-13]。地理加权回归模型(geographically weighted regression,GWR)为局部模型,在不同地理空间位置上的估计参数不同,是一种非常有效揭示空间非平稳性和空间依赖的方法[14]。因此,采用GWR探究影响植被变化的地理因子及其空间分布的异质性,对于揭示植被变化的机制与机理,指导不同区域因地制宜的生态恢复与重建具有重要的指导意义。
宁陕蒙接壤地区的宁夏宁东能源化工基地、内蒙古鄂尔多斯市和陕西榆林市在地理上构成了一个三角地带[15],统称能源“金三角”地区。该区域地理位置优越,是国家特大级能源基地和全国煤炭外输基地,2013年该地煤炭外运量占煤炭总产量的80%,在我国能源发展中发挥举足轻重的作用。随着“一带一路”倡议的提出,中国与欧盟、中亚地区能源合作的展开,能源“金三角”地区将成为重要的能源通道。目前,国内学者对于能源“金三角”地区的研究较多,主要涉及发展战略、空气质量、水资源保护和植被动态监测等方面[15-18]。但是关于定量研究该地区植被变化特征及其驱动因素的报道较少,特别是应用GWR分析不同地理位置上不同因子对植被变化影响的研究更少。在能源基地建设过程中，生态环境保护尤为重要,能源“金三角”地区位于相对封闭的干旱半干旱内陆,生态系统状态及其变化对于中国西部地区环境有重要影响。据此,以能源“金三角”地区作为研究区域,基于MODIS-NDVI数据,对2005—2015年的植被动态变化趋势、空间分布等进行分析,基于GWR模型对植被动态变化的影响因素进行研究,以期为该地区能源战略规划和生态环境保护提供科学依据,从而有助于实现产业与生态的同步发展。
1 研究区概况
能源“金三角”地区位于36°24′～ 40°51′ N,104°17′～111°27′ E之间(图1)，地处国家中、西部过渡地区的鄂尔多斯盆地,地广人稀,总面积13.35万km2,人口约553.8万人。
图1 研究区概况
该区属于内陆半干旱气候区,多年平均降水量在400 mm以下,年平均蒸发量在2 000 mm以上;水资源贫乏,植被覆盖率低,水土流失严重;西部沙丘沙地广布,东部沟壑梁峁发育,是我国典型的生态环境脆弱区。
2 数据来源与研究方法
2.1 数据来源与处理
2.1.1NDVI数据获取与处理
NDVI数据来源于美国国家航空航天局(NASA)的MOD13Q1 NDVI数据产品,空间分辨率250 m,时间分辨率16 d,时间跨度为2005—2015年,考虑能源“金三角”地区植被物候特征,下载每年植被生长旺盛期(7—9月)影像共176景。首先利用MRT(MODIS reprojection tools)完成图像的空间拼接;然后,为了消除异常值的影响,利用ENVI 4.8软件采用最大合成法(maximum value composites,MVC)[19]生成年NDVI数据。年最大NDVI能够有效反映区域植被的覆盖状况,将研究区域各像元年最大NDVI的平均值作为当年的NDVI。2013年NDVI均值达最大值,因此,以2013年作为植被变化节点,用2个时段NDVI的差值来表征研究区植被变化情况,2005—2012年为植被增长阶段,2013—2015年为植被减少阶段。
2.1.2影响因素指标获取与处理
利用ArcGIS软件的空间统计工具探测NDVI变量的空间相关性。NDVI指数时空格局受自然因素和社会经济因素相互影响,根据文献[20-22],综合考虑自然环境与社会经济等因子的基础上,兼顾数据的可获取性,结合研究区域特点,并满足空间化、定量化等要求,最终选取高程、坡度、土壤黏粒含量、多年平均气温、多年平均降水、距煤矿区距离及距道路距离共7个影响因素。其中,DEM数据来自于中国科学院资源环境科学数据中心中国海拔高度DEM(SRTM 90 m)空间分布数据,并据此生成高程和坡度数据;土壤数据下载自中国科学院资源环境科学数据中心;气象数据采用中国气象科学数据共享服务网提供的中国地面气候资料月值数据集的月平均温度和月降水量资料,利用ArcGIS软件对气象数据进行Kriging空间插值;道路数据通过数字化《中华人民共和国分省系列地图》得到〔GS(2010)1541号,中国地图出版社〕。各影响因子统一为5 km×5 km规则网格,气象、土壤数据与DEM等栅格数据通过ArcGIS软件空间分析模块中的区域分析方法将其单元统一;距离因子通过 ArcGIS软件中的邻域分析统计网格平均值;采用离差标准化方法对上述数据进行标准化处理。
2.1.3多重共线性检验
回归分析时,变量间相关性程度高、信息重叠,容易产生多重共线性问题,从而影响回归模型的正确估计。方差膨胀因子(VIF)等检验方法可判断多重共线性,剔除冗余变量。研究利用EVIEWS软件中VIF指标对影响因子进行多重共线性进行检验,排除潜在多重共线性问题。一般认为,VIF≥10,存在强共线性;5≤VIF<10,存在中度共线性;2≤VIF<5,存在轻度共线性[23]。
2.2 研究方法
2.2.1NDVI趋势分析
文章采用一元线性回归方程的斜率来模拟研究区NDVI的变化趋势,即对随时间变化的NDVI进行回归分析，生成年平均NDVI数据,该图像中的每个像元都有11 a的相应时间序列数值,对这些数值进行回归分析,来揭示每个栅格NDVI在11 a时空序列中的演变趋势,并得到NDVI演化趋势图。计算公式为
(1)
式(1)中,变量i为年份序号,取值为1～11;INDV，i为第i年的NDVI值;n为研究时段的年份数，a;s为2005—2015年NDVI的线性回归斜率。若s>0,说明此像元NDVI呈增加趋势,反之则减少。为了解变化趋势的显著性,对变化趋势进行F检验。结合s值和F检验结果,将NDVI年际变化趋势分为5个等级[24]:显著减少(s<0,P≤0.01);不显著减少(s<0, 0.01
2.2.2GWR模型
GWR模型是对传统回归的拓展,其权为观测值所在的地理空间位置与回归点的地理空间位置之间的距离函数,模型表达式为
(2)
式(2)中,yi为样本i的NDVI变化拟合值;(ui,vi)为第i个样本空间单元的地理中心坐标;β0(ui,vi)为第i个样本的常数项估计值;βj(ui,vi)为第i个样本的第j个回归参数,是关于地理位置的函数;xij为第j个自变量在样本i的值;εi为服从均值为0的独立正态分布的误差。
由于GWR模型在诊断方面效果较弱,因此在进行地理加权回归分析之前需要先进行普通最小二乘法(ordinary least squares，OLS)线性回归以确保模型的准确性[25]。普通最小二乘法为全局性线性回归模型,即在各区域用全部自变量估计因变量的值,基本原则是最优拟合直线,使各点到直线距离的平方和最小,回归方程表达式为
(3)
式(3)中,yi为因变量,是自变量xij(j=1, 2,…,k)的一个线性组合,即NDVI变化的拟合值;β0为常数项;βj为第j个自变量的回归系数;i为观测值个数;εi为符合正态分布的随机误差项。
研究采用Gauss函数,利用GWR 4.0工具实现建模,根据赤池信息量准则(akaike information criterion,AIC，CAI)、AIC值最小原则确定最优带宽。将研究区域划分为5 km×5 km的栅格,将2005—2012年7—9月平均NDVI差值与2013—2015年7—9月平均NDVI差值作为GWR模型的因变量,7个影响因素作为自变量。将每个栅格回归系数绝对值最大值所对应的影响因素作为网格植被变化的主要影响因素,进而分析主导因素的空间异质性。
2.2.3模型评价
为比较OLS模型和GWR模型拟合效果,选用赤池信息量准则检验对比模型显著性,表达式为
(4)
此外,模型中残差平方和(residual sum of squares)为评价回归模型的拟合效果指标,残差平方和越小,拟合优度越高。利用Moran′sI对GWR模型误差项进行空间自相关分析,Moran′sI值接近0表示模型的残差随机分布,拟合效果较好;反之,模型拟合效果较差。
3 结果与分析
3.1 NDVI时空变化分析
3.1.1NDVI时间变化特征
由图2可见，11 a间研究区和煤矿区NDVI均值虽有波动,但整体呈现不断增加趋势,增长率分别为0.083·(10 a)-1和0.090·(10 a)-1(P<0.05),表明植被处于改善状态。11 a间,年均NDVI变化经历了2006—2007和2011—2012年2次较大波动,增长率分别达0.320·(10 a)-1和0.335·(10 a)-1。从煤矿区内外来看,研究区域各年NDVI均值高于煤矿区内NDVI均值,说明煤炭开采对植被造成一定不良影响。研究区域NDVI均值波动幅度在0.337～0.458间;煤矿区内NDVI均值波动幅度在0.242～0.383间;NDVI数值相对较低,波动范围较小。2006与2007年煤炭产量平稳,煤矿区内植被长势良好;2008—2012年煤炭产量直线上升,研究区与煤矿区植被覆盖基本稳定。值得注意的是,2011—2012年煤炭继续增产的同时,研究区与煤矿区NDVI增加,到2013年达最大值0.458;2013—2015年煤炭产量稳定,植被覆盖却急剧下降,说明植被生长不仅受资源开发影响,仍有其他影响因素需要进一步探索。据此可以把研究区植被变化划分为振荡上升期(2005—2012年)和下降期(2013—2015年)2个阶段,来研究植被变化的主导因素。
图2 2005—2015年NDVI值变化
3.1.2NDVI空间分布特征
利用ArcGIS软件将研究区2005—2015年最大化合成的NDVI影像进行逐像元平均,得到研究区11 a平均NDVI空间分布图(图3),并按照自然分段法分5级显示。能源“金三角”地区NDVI东高西低,呈现由东南向西北逐渐递减的空间分布特征。NDVI最小的区域主要分布在鄂尔多斯市杭锦旗和乌审旗,为沙漠与沙地覆被类型,占总面积的17.21%;植被覆盖较好的区域主要分布在榆林市南部黄土高原丘陵沟壑区,多为森林草原覆盖,NDVI值为0.45～0.58,占总面积的34.49%;NDVI值大于0.58的区域仅占3.86%。
图3 2005—2015年多年平均NDVI空间分布Fig.3 Spatial distribution of average NDVI during the period from 2005 to 2015
2005—2015年能源“金三角”地区NDVI年际变化趋势的显著性分析(图4)表明,显著减少、不显著减少、稳定、不显著增加和显著增加的面积分别占研究区总面积的0.27%、0.40%、72.23%、13.43%和13.68%。NDVI变化呈增加趋势的区域(27.11%)远大于减少区域(0.64%),呈稳定趋势面积比例高达72.23%,表明11 a来研究区NDVI整体呈现稳中增加的趋势。
图4 2005—2015年NDVI变化趋势显著性空间分布
从显著性空间分布来看,NDVI减少(显著和不显著)区域在研究区内零星分布,增加(显著和不显著)区域主要分布在杭锦旗北部,榆林市东部神木县、榆阳区、横山县、佳县、米脂县、吴堡县、清涧县等,以及宁东的灵武市和红寺堡开发区。杭锦旗北部为黄河生态旅游区,植被增加显著;榆林市植被生长较好的原因可能是陕西中北部退耕还林和防风固沙政策实施较早[26],政策的推行有利于植被恢复与生长。而未发生显著变化的区域植被覆盖相对较少,库布奇沙漠、毛乌素沙地等自身环境相对不利于植被生长,生态工程效果还没有得到很好的体现,是未来工作的重点和难点。
3.2 NDVI变化的驱动因素分析
3.2.1OLS模型及结果
以NDVI增长值作为因变量,与自变量高程、坡度、土壤黏粒含量、多年平均气温、多年平均降水、距煤矿区距离及距道路距离7个影响因素进行OLS回归,结果如表1所示,各因素VIF值都不超过2,表明自变量中仅存轻度共线性,但2013—2015年距煤矿区距离因素未通过10%水平下的显著性检验,因而将其余6个因素纳入OLS模型。
比较模型中自变量的标准化系数可以发现,2005—2012年对植被增加影响程度最大的是气温,其次是土壤黏粒含量与坡度(表1)。OLS模型结果表示,其他条件不变的情况下,气温每增加1个标准化单位,NDVI增加0.060个单位;土壤黏粒含量每增加1个标准化单位,NDVI增加0.003个单位;2013—2015年,对植被增加影响程度最大的是坡度,坡度每增加1个标准化单位,NDVI下降0.023个单位。OLS模型总体上说明地形气候等自然要素对NDVI增长值的影响程度大于人类社会要素。
3.2.2GWR模型及结果
应用GWR模型分析的前提条件是变量间的空间相关性确实存在,因此,在建模和分析之前使用Moran′sI值对NDVI变量的空间相关性进行预检验。利用ArcGIS软件的空间统计工具开展计算,结果显示研究区域2005—2012与2013—2015年Moran′sI值分别为0.825和0.851,揭示该区域NDVI具有较强的集聚模式特征。
在2个时段中,GWR模型具有最大R2、最小AIC 和残差平方和,表明基于GWR模型的NDVI变化驱动因素的拟合效果优于OLS模型。如表2所示,对比2个模型的拟合效果,GWR 模型回归的拟合优度达0.74和0.77,调整的拟合优度为0.72和0.76,而OLS模型回归的拟合优度仅为0.23和0.29。这意味着GWR模型解释了NDVI变化的74%和77%,OLS模型最多只能解释NDVI变化的29%。GWR模型可在局部水平上对NDVI变化量进行解释,其解释程度较OLS模型提高46百分点以上。
表1OLS模型参数估计与检验结果
Table1ParameterestimationandtestresultsintheOLSmodel
时段因素参数估计值标准误差t值P值方差膨胀因子VIF 2005—2012年常数项-0.405 000 0000.018-22.1430.000高程-0.000 014 1000.000-3.7020.0001.126坡度0.002 000 0000.0012.1470.0321.482土壤黏粒含量0.003 000 0000.00016.7650.0001.108多年平均气温0.060 000 0000.00231.4530.0001.290多年平均降水-0.000 119 0000.000-9.1390.0001.642距煤矿区距离-0.000 000 1190.000-5.3790.0001.480距道路距离0.000 000 4900.0007.5500.0001.085 2013—2015年常数项0.251 000 0000.01320.0280.000高程-0.000 072 2000.000-20.8660.0001.099坡度-0.023 000 0000.001-29.1760.0001.222土壤黏粒含量-0.002 000 0000.000-16.0780.0001.084多年平均气温-0.021 000 0000.001-15.2930.0001.894多年平均降水0.000 015 3000.00012.0730.0001.876距道路距离0.000 000 1380.0002.3060.0211.092
表2OLS模型与GWR模型的拟合效果比较
Table2IndicesforcomparingOLSandGWR
时段模型带宽残差平方和AIC值拟合优度R2调整拟合优度R2 2005—2012年GWR300.006.51-25 637.420.740.72OLS19.28-19 238.910.230.23 2013—2015年GWR298.995.11-27 273.980.770.76OLS15.96-20 463.580.290.29
对所有网格的标准化残差进行空间自相关性检查,得到空间自相关系数分别为0.003和0.006,Z值统计量分别为0.649和1.083,说明地理加权回归预测值标准差在空间上随机分布,模型整体模拟效果较好。因此,利用GWR模型对能源“金三角”地区植被变化的影响因素进行回归分析合理可行。
与OLS模型是一个全局或平均意义的估计值不同，GWR模型的计算结果中,每个格网都有特定的回归系数。统计各影响因素的回归系数,得到平均值、最小值、最大值、中位值和正负影响所占比例(表3)。统计结果表明,各影响因子对格网内的NDVI变化既存在正向影响又存在负向影响,2005—2012年,正值比例和负值比例最高的分别是坡度(73.21%)和距矿区距离(64.12%);2013—2015年正值比例和负值比例最高的分别是多年平均气温(57.81%)和坡度(79.82%)。为进一步验证各影响因素回归系数是否具有空间非平稳性,对回归系数进行空间自相关性分析,得到表3中各因素回归系数的莫兰指数与Z统计量,说明在1%显著水平下,各项系数具有显著的正向空间自相关。因此,GWR模型可以较详细地分析不同影响因素在不同地理位置上对植被变化的定量影响。
3.2.3植被变化的主导因素空间格局分析
提取每个格网对植被生长影响最大的因素(图5),分析造成NDVI变化差异的原因,发现影响植被变化的主导因素随着空间地理位置不同而变化,且各个影响因子与植被变化既有正相关又有负相关。研究表明,多年平均气温与降水是研究区大部分地区的主导因素,在2个时段分别占总面积的86.29%和88.68%;但在空间分布上有所差异,2005—2012年研究区北部主要受多年平均气温影响,中部与南部主要受多年平均降水影响;2013—2015年研究区西部受多年平均气温影响,东部依然主要受多年平均降水影响。分析认为,能源“金三角”地处黄土高原地区,该地区十分缺水,土壤保水性较差,侵蚀严重,土壤含水量低且主要分布在深层土壤[27],植被生长所需水分主要来自降水;因此,降水给该地植被带来较明显的直接影响。研究区降水量由东南向西北逐渐减少,西部多沙丘沙地,东部沟壑梁峁发育,丘陵沟壑区域内部由于环境小气候的影响,植被变化的影响因素以降水量为主[28]。研究区日照充足,西北部年辐射总量较大,光热资源丰富，有利于植被生长[29],但毛乌素沙地气温升高使蒸散增加，导致植被发生干旱胁迫。在平原与丘陵地貌相间、地形较复杂的地区，地形因素在2013—2015年较2005—2012年表现出更显著的影响。2个时段中土壤黏粒含量没有表现出对植被变化的主导作用,已有研究显示,地形因子对土壤含水量、颗粒状况等有重要影响[30],因此,地形因素与土壤因素同时存在时,地形因素成为影响植被变化的较主要因素。研究中采用的交通因素是密切影响区域经济发展情况、劳动力与资本流通情况的铁路、省道、国道及高速公路等道路,该因素的主导作用仅占研究区总面积的1%。
2005—2012年能源“金三角”地区能源开发量快速增加,煤炭产量逐年增长,煤矿区因素在神东矿区、榆横矿区、榆神矿区、灵武矿区、马家滩、积家井等矿区均表现出对植被变化的主导作用;2013—2015年,该地区实施煤炭资源整合,煤炭产量稳定,对植被变化影响不显著。由此,影响研究区植被变化的主导因素在时空上存在较显著的差异性。
表3GWR模型回归系数统计
Table3StatisticsofregressioncoefficientsintheGWRmodel
时段因素 平均值最小值最大值中位值正值比例/%负值比例/%Moran′s I值Z值 2005—2012年常数项0.16-1.302.680.1178.0121.99高程0.00-0.120.140.0051.1348.870.98110.70坡度0.01-0.050.130.0173.2126.790.97110.26土壤黏粒含量0.00-0.040.040.0045.9654.040.97109.72多年平均气温-0.01-0.941.610.0153.7046.300.96108.56多年平均降水-0.02-0.521.02-0.0341.7058.300.96108.55距煤矿区距离-0.02-0.290.37-0.0235.8864.120.97109.46 距道路距离0.00-0.050.040.0056.8343.170.97108.77 2013—2015年常数项-0.02-0.821.10-0.0243.4856.52高程-0.02-0.130.11-0.0130.8069.200.98110.96坡度-0.01-0.080.05-0.0120.1879.820.98110.40土壤黏粒含量0.00-0.050.070.0031.0968.910.97108.97多年平均气温0.05-0.841.020.0257.8142.190.99111.15多年平均降水-0.03-1.791.26-0.0146.3853.620.98110.28距道路距离0.00-0.060.050.0047.2652.740.98110.77
图5 2005—2015年植被变化主导因素的区域差异分布Fig.5 Regional differences of dominant factors on vegetation dynamics
4 结论与建议
基于MODIS传感器的MOD13Q1数据集探讨了能源“金三角”地区2005—2015年间NDVI的时空变化特征,定量化分析了影响NDVI变化的主要驱动因素及其空间非平稳性。研究表明:
(1)时间上,2005—2015年研究区平均NDVI整体呈现波动增加趋势,增长率为0.083·(10 a)-1(P<0.05);煤矿区植被生长水平低于全区生长水平,NDVI增长率为0.090·(10 a)-1(P<0.05)。空间上,NDVI呈现由东南向西北递减的分布格局;趋势上,NDVI变化呈增加趋势的区域(27.11%)远大于减少区域(0.64%),显著增加区域主要分布在榆林市东部。
(2)NDVI及其影响因子属于地理空间数据,具有空间自相关性和空间非平稳等特征,采用GWR模型对能源“金三角”地区植被变化的影响因素进行回归分析合理可行,且提高了46百分点以上的可靠性。以空间同质性为前提假设的全局模型OLS不能很好地解释NDVI与影响因素的空间非平稳性,OLS模型虽有一定的解释能力,但由于没有考虑到空间自相关的影响，仍不能全面地解释NDVI与影响因素的空间分布。与2005—2012年相比,2013—2015年全区NDVI显著减少,且具有更强的空间聚集性(Moran′sI值为0.851),但变化程度存在区域差异。
(3)2005—2015年能源“金三角”地区植被变化整体受气候等自然因素的影响较大,影响植被变化的主要因素在时间与空间上具有差异性。2005—2012年研究区北部的主要影响因素为多年平均气温,中部与南部的主要影响因素为多年平均降水;2013—2015年研究区西部受多年平均气温影响,东部主要受多年平均降水影响。
人类活动对植被变化具有双重作用。一方面,煤矿区生态环境的重建、退耕还林还草、土地荒漠化防治、禁牧、种草种树的力度对植被的长势起到促进作用;道路景观与绿化带的建设也对植被保护起到促进作用。另一方面,能源“金三角”地区属于典型的环境脆弱区,某些不可逆的活动使得煤矿采空区植被依然出现退化现象;经济发展伴随的道路工程、道路网络对植被格局呈干扰作用。然而,由于不同地理位置自然条件、开发建设程度、保育成效等方面的差异,植被对其响应也具有空间差异性。
将地理加权回归模型应用于植被变化与影响因素研究,判定了不同空间位置模型的拟合度,得出各影响因子随空间位置变化的权重;与传统的植被覆盖等研究相比,更能整体、直观地揭示该地区内植被的空间特征。基于以上结论,不同区域影响植被变化的主导因素存在差异,建议:(1)可采用分区治理的方式,综合考虑植被变化主导因素的区域差异,针对不同情况制定差异化植被恢复与保护方案;(2)将资源开发与环境优化相结合,构建有利于生产发展和生态保护相协调的体制机制。在下一步工作中,该地区的开发应协调好与脆弱的生态环境之间的关系,保障能源资源可持续利用与区域经济社会环境长期平稳健康发展。
