植被是陆地生态系统的主题,是保护生态系统的重要屏障。2014年《河南省主体功能区规划》确定构建全省以“四区三带”为主体的生态安全战略格局,其中伏牛山生态区为河南省重点生态功能区。因此,本文选择伏牛山生态区作为本文的研究区。研究区位于河南省西南部,见图1。西与陕西省接壤、东至鲁山县、南召县,南至淅川、内乡县与南阳盆地的北缘相连、北至陕县、宜阳县,面积约为34 471.5 km²。研究区温凉湿润,少旱区,地处中国南北气候的过渡带,属于北亚热带向暖温带的过渡地带。伏牛山山区气温:春季,4月份平均气温在5℃以上,平均每月增温为5～6℃,是全年增温最快的季节; 夏季,气温达到最高,最热月7月平均气温为26.5～28.5℃; 秋季,9月基本上与4月相似,是气温日较差最大的季节。伏牛山山区年平均气温在12℃左右,山区年平均降水量在800～1 100 mm,7月、8月降水量为全年最多。植被属暖温带落叶阔叶林向亚热带常绿阔叶林的过渡类型,生物多样性十分丰富。主要植被类型分为7种,分别是针叶林,阔叶林,针阔叶混交林,竹林,灌丛及灌草丛,草甸,沼泽及水生植被^[17]。伏牛山地区森林植被保存完好,森林覆盖率达42%。

图1 研究区分布概况

植被覆盖度的变化受气候影响较大。北半球植被覆盖度与20世纪80,90年代相比,增长趋势变缓,气温与降水等气象因子是影响植被覆盖度变化的主要驱动力^[18-21]。

本研究采用的气象基础数据来源于河南省气象局和中国气象科学数据共享服务网^[22]。气象基础数据包括了研究区域2000—2018年内共14个气象站的4-9月的月均温和月降水量数据。在ArcGIS中进行矢量化,并保证其坐标系与遥感数据的坐标系一致。运用克里金方法对气温和降水量数据进行克里金插值可以得到研究区域栅格化的月均温和月降水量的分布。研究区气象数据统计如图2。研究区的降水与高温是同步的,属于雨热同期的气候,研究区平均气温4月最低,为16℃,7月最高为26.4℃,而研究区的月降水量4月最低为20.3 mm,7月最高为217 mm。

图2 研究区2000-2018年月降水量和月均温变化

选用2017年7月的landsat8影像作为数据精度验证,首先进行辐射定标,然后基于FLAASH大气校正模块对影像进行大气校正,将landsat8影像变换坐标系、重采样为500 m,并进行植被覆盖度的计算。随机选取54个点,分别提取landsat8的植被覆盖度与同期的植被覆盖度进行线性拟合,得到的拟合方程为y =0.496 x+0.380 1,R²为0.71。说明利用MODIS数据基于像元二分模型计算的植被覆盖度可信度较高。

利用基于标准化植被指数(NDVI)^[23]的像元二分模型对伏牛山研究区植被覆盖度(VFC)进行估算:

NDVI=NDVI_veg+NDVI_soil(1)

式中:NDVI_soil为完全是裸土或无植被覆盖区域的NDVI值; NDVI_veg表示完全被植被覆盖的像元的NDVI值。模型中仅把像元划分为两类:裸土和纯植被,这样可以削弱大气,植被类型等因素的影响。

NDVI=VFC×NDVI_veg+(1-VFC)×NDVI_soil(2)

按照所占比例分析计算,进行等式变化可得:

VFC=((NDVI-NDVI_soil))/((NDVI_veg-NDVI_soil))(3)

由于图像中不可避免的存在噪声,因此NDVI的极值并不一定是NDVI_max与NDVI_min。在取值时需要选取一定置信区间内的最大值与最小值。选取累积概率为1%,99%的NDVI值作为NDVI_min和NDVI_max。基于像元二分模型对研究区植被覆盖度进行估算,并将植被覆盖度进行等级划分,分为5个等级(表1)。

表1 植被覆盖度等级划分

对伏牛山地区的植被覆盖度变化情况采用马尔科夫模型进行分析,得到一个转移矩阵,可以定量的表明不同等级的植被覆盖度之间的变化情况。

V_ij=[V₁₁ V₁₂ … V_1n

V₂₁ V₂₂ … V_2n

   

V_n1 V_n2 … V_nn](4)

V_ij表示研究区第i种等级的植被覆盖度和第j种等级的植被覆盖度之间转化的面积。n为覆盖度的等级数量。i和j的取值为1,2,…,n。

利用逐像元拟合年均VFC(植被覆盖度)的变化斜率得到植被覆盖度的时空变化数据,通过统计时空变化数据可以对每年的植被覆盖度做出在年序上的时间变化趋势分析。植被覆盖的变化趋势可以综合反映研究区内植被的时空发展演变特征^[24]。利用一元线性回归分析法逐像元计算2000—2018年的趋势斜率:

Slope=(n∑ⁿ_i=1(iVFC_i)-∑ⁿ_i=1i∑ⁿ_i=1VFC_i)/(n∑ⁿ_i=1i²-(∑ⁿ_i=1i)²)(5)

式中:n是时间序列的长度(n=19); VFC_i是第i年的VFC值; Slope代表研究时段内每个像元的植被覆盖度的变化趋势,当Slope>0时,代表该像元处19 a间植被覆盖度呈增长趋势,反之则表示减少。显著性趋势检验采用F检验,若相关系数置信度达到95%,99%以上即p<0.05,p<0.01,则认为变化趋势分别达到了显著和极显著水平。根据检验结果将变化趋势分为5个等级:极显著增加(p<0.01,Slope>0),极显著减少(p<0.01,Slope <0),显著增加(0.01<p<0.05,Slope>0),显著减少(0.01<p<0.05,Slope<0),不显著(p>0.05)。

自然界中具有长程依赖性的时间序列是普遍存在的,并在水温,地球化学,气候,地质和地震等领域广泛运用^[25]。而Hurst指数是描述长程依赖性的一种方法。Hurst指数可以反映出序列自相似性以及序列发展的相关强度^[26]。基于重标极差(R/S)分析法计算得到的Hurst指数可以反映序列自相似性及序列发展的相关强度,并能够较好的判别时间序列数据的自相关性。若0<H<0.5,则表明该时间序列具有长程依赖性,即未来变化趋势与过去的趋势相反,称其为反持续性。H值越趋近于0,其反持续性越强。若H =0.5,则时间序列是相互独立的,说明时间序列不存在长期相关的特性。若0.5<H<1,说明时间序列是一个持续性序列,未来变化趋势与过去的趋势相同即同向持续性。H值越趋近于1,说明其同向持续性越强。

偏相关分析是反映多个变量间的相关程度,将伏牛山地区的植被覆盖度分别与气温和降水量逐像元的进行偏相关分析。偏相关系数的计算公式如下:

式中:r_xy1,y2表示将变量y₂固定,对变量x和变量y₁进行偏相关分析得到的偏相关系数; r_xy1是变量x和变量y₁的简单相关系数; r_xy2是变量x和变量y₂的简单相关系数; r_y1,y2是变量y₁和变量y₂的简单相关系数。r_xy1,y1取值范围是[-1,1],当r_xy1,y1>0则表示两个变量间呈正相关,r_xy1,y2<0则表示两个变量间呈负相关。r_xy1,y2越接近-1或1,说明两个变量间的偏相关性越高。

选取2000年、2010年、2018年为代表年,对其进行分析统计,可以得到研究区植被覆盖度等级分布图(图3),卢氏、栾川、嵩县、西峡、南召、内乡地区的较高覆盖度和高覆盖度的像元的个数较多,即植被覆盖度高于0.6的地区主要集中在这些地区,说明这些地区植被覆盖水平较高,同时也反映了伏牛山山区植被覆盖度较高。对于伏牛山周边的地区,即灵宝,洛宁,宜阳,汝阳等地区,其靠近伏牛山一侧的植被覆盖度要高于远离伏牛山一侧的植被覆盖度。即研究区的植被覆盖度总体上呈现中间高,四周低的特点,与伏牛山在研究区内的分布情况相一致,说明人类的生产生活等行为会对植被覆盖度的变化产生较大的影响。

图3 2000、2010、2018年研究区植被覆盖度等级分布

根据同一地区不同时相的土地覆盖现状可以求得一个转移矩阵。分析转移矩阵,能够得到两个时相,不同土地覆盖类型之间相互转化的情况^[27]。对研究区不同等级的植被覆盖度计算其转移矩阵并分析,可以得到植被覆盖度的转移矩阵见表2—3。

由表2可知2000—2010年研究区的植被覆盖度有59.4%的区域没有发生变化。低覆盖度区转出面积为116.25 km²,转入面积为95.25 km²转出面积大于转入面积,表明低覆盖度区域在减少,当地植被覆盖度水平在上升。较低覆盖度区转出面积为2 902.25 km²,转入面积为273.75 km²,其中转出到中覆盖度区的面积占较低覆盖度区转出面积的92.25%,表明较低覆盖度区变化较大,变化趋势是向良性发展的。中覆盖度区转出面积为6 539.5 km²,转入面积为2 831.5 km²,转出到较高覆盖度区的面积占中覆盖度区转出面积的97.3%,表明研究区内生态环境较好,植被覆盖度呈增长的趋势。较高覆盖度区转出面积为158.5 km²,转入面积为7 650.5 km²。转入面积大于转出面积,较高覆盖度区的面积有大幅度的增加。表明当地生态环境得到改善,大部分劣等覆盖度区转为较高覆盖度区。高覆盖度区转出面积为1 147 km²,转入面积为12.5 km²。高覆盖度区的转出面积大部分转为了较高覆盖度区。表明较高覆盖度区和高覆盖度区相互转化较多。

表2 研究区2000年、2010年植被覆盖度转移矩阵

表3 研究区2010年、2018年植被覆盖度转移矩阵

由表3可知2010—2018年研究区的植被覆盖度有66.32%的区域没有发生变化。低覆盖度区转出面积为118 km²,转入面积为271.5 km²,转出面积小于转入面积,低覆盖度区域的面积有所增加。较低覆盖度区转出面积为393.75 km²,转入面积为926.25 km²,其中由中覆盖度区转为较低覆盖度区占较低覆盖度区转出面积的85.4%,劣等植被面积有所增加。中覆盖度区转出面积为3 202.25 km²,转入面积为2 215.5 km²,中覆盖度区面积有所下降。较高覆盖度区转出面积为7 823 km²,转入面积为2 439 km²,转出面积大于转入面积,且转出面积是转出面积的3.2倍,表明较高覆盖度区面积下降严重。高覆盖度区转出面积为102.75 km²,转入面积为5 787.5 km²,转入面积为转出面积的56倍,表明高植被覆盖度区面积增加明显,生态环境日益改善,生态治理有一定的效果。

综上可知,2000—2010年低覆盖度区、较低覆盖度区及中覆盖度区的面积在减少,较高覆盖度区和高覆盖度区的面积在增加,研究区的植被覆盖度呈上升趋势。2010—2018年较低覆盖度区和高覆盖度区的面积增加显著,而较高覆盖度区的面积有所下降,低覆盖度区和中覆盖度区的面积有所变化,但变化不大,研究区的植被覆盖度整体呈上升趋势。即近19 a研究区内的较多低水平植被覆盖度区转化为高水平植被覆盖度区,说明了当地十分重视保护生态环境。

由像元二分模型计算得到研究区2000—2018年植被覆盖度分布图。使用一元线性回归分析算法逐像元的分析覆盖度变化并进行显著性检验,得到研究区覆盖度变化趋势图,见图4。研究区植被覆盖度整体呈逐年上升的趋势,呈增加趋势的地区占整个研究区的68%,5%的区域其植被覆盖度呈不同程度的下降。对Slop进行F检验,并对研究区的植被覆盖度不同变化趋势进行等级划分,统计不同等级的植被覆盖度区域所占比重和其面积。73%的研究区域通过了显著性检验,其中极显著增加的区域和显著增加的区域所占比重分别为47%和21%,通过显著性检验且呈下降趋势的区域(极显著减少区域占比2%,显著减少3%)。下降的地区主要集中在灵宝市城区,洛宁东部地区,宜阳东部地区,栾川县城区,内乡县城区以及镇平县城区,说明这些地区19 a来城镇化水平不断提升。未来为了更好的保护研究区内的森林等自然资源,需要控制城镇用地的扩张。

基于重标极差(R/S)分析法逐像元计算研究区的Hurst指数,并绘制研究区植被覆盖度的Hurst指数空间分布图。研究区的Hurst指数范围为0～0.98,平均值为0.47,其中Hurst指数大于0.5和小于0.5的像元分别占研究区的65%和35%。将Hurst指数与通过显著性检验的Slop结合,可以对研究区内植被覆盖度的未来趋势进行预测,预测结果如图5所示。研究区未来覆盖度增加的地区仅占整体研究区的26.5%。未来覆盖度减少的区域占比为46.3%,接近一半的区域其植被覆盖度将延续减少的发展趋势。

图4 2000-2018年研究区植被覆盖度变化趋势显著性

图5 研究区植被覆盖度未来变化趋势

统计各变化趋势所占面积和比重可以发现,未来极显著增加地区的面积为6 103.5 km²,所占比重为17.1%,主要分布在栾川县,汝阳县,卢氏县,西峡县和南召县的伏牛山山区,由于对伏牛山山区的森林进行了多年的保护,这些地区的植被覆盖度未来呈极显著增长的趋势。对于未来呈显著增加趋势的地区,其主要分布在伏牛山山区,所占面积为3 591.5 km²,占比仅为9.4%,这些地区植被覆盖度未来转向良性发展或者继续保持当前良好的发展状态。而未来极显著减少和未来显著减少的地区分布范围较广,其中未来极显著减少占比高达32.4%,表明未来覆盖度呈下降趋势的地区占比较大。面对脆弱的生态环境,需要坚持退耕还林的政策,应该重视未来呈极显著减少和显著减少地区的生态环境治理。

由于气温和降水都会对研究区内的植被覆盖度的变化产生影响,将研究区内的植被覆盖度分别与月降水量和月均温进行偏相关分析,先选取代表年进分析,2000年植被覆盖度与降水和均温的平均偏相关系数分别为0.24,0.18,呈正相关的地区分别为78%和64%。2010年植被覆盖度与降水和均温的平均偏相关系数分别-0.13,0.54,呈正相关的地区分别为36%和84%。2018年植被覆盖度与降水和均温的平均偏相关系数分别为0.14,0.38,呈正相关的地区分别为60%和75%。然后对2000—2018年共19 a间的植被覆盖度与降水和均温进行偏相关分析,见图6,植被覆盖度与同期降水量和气温的平均偏相关系数分别0.067,0.097,呈正相关的所占比重分别为56.18%和63.97%。研究区内的植被覆盖度与气温的相关程度更紧密。进一步对研究区内的森林区域进行研究,发现研究区内的森林区域的植被覆盖度与气温的偏相关系数为0.12,呈正相关的区域占森林区的70%,与降雨的偏相关系数为-0.13,呈正相关的区域占森林区的27%,表明气温是影响森林覆盖度的主要因子(图7)。

图6 2000-2018年研究区植被覆盖度与降水量和温度的偏相关系数

图7 2000-2018年研究区森林地区植被覆盖度与降水量和温度的偏相关数