基于像元二分模型的伏牛山地区植被覆盖度变化

张成才, 娄 洋, 李 颖, 姬兴杰, 董萌佳

(1.郑州大学 水利科学与工程学院, 郑州 450001; 2.中国气象局 河南省农业气象保障与应用技术重点开放实验室, 郑州 450003; 3.河南省气象科学研究所,郑州 450003; 4.河南省气候中心, 郑州 450001; 5.河南省白沙水库管理局, 河南 禹州 461670)

植被作为地表生态系统重要的成分之一，与土壤、水、大气等构成一个生态循环系统[1]。植被覆盖受气候变化和人类活动的影响较大，开展植被覆盖研究揭示生态系统的变化特征成为了一个重要的研究领域。植被覆盖度是指某一区域内植被(包括叶，茎等)投影面积与该地域面积之比，它可以清晰的反映研究区的植被情况，是一项重要的生态环境参数[2]。越来越多的学者选择植被覆盖度作为量化指标进行植被时空变化和区域生态的研究。植被覆盖度反演主要有地面实测和遥感估算两种方法。地面实测目前所采用的方法有目测法，采样法，仪器法等，常用于小尺度测量。由于其精度较高，一般作为检验遥感估算覆盖度的手段。遥感估算植被覆盖度逐渐成为分析植被覆盖度变化的主要方法[3]。它适用于大尺度反演植被覆盖度，并且具有良好的时间和空间连续性，可以从多维度研究分析植被覆盖度的差异及其变化的规律[4]。遥感估算目前所采用的方法有回归模型法、植被指数法、混合像元法以及机器学习法等[5]。

目前，植被覆盖度的研究成果主要包括了对植被覆盖度动态变化的研究[6-7]、对植被覆盖度变化与生态环境因子的关系进行探究[8-9]、以及对植被覆盖度监测的研究[10-11]。赵舒怡等[12]对华北地区植被覆盖度与干旱条件进行了相关分析。张晓东等[13]利用TM卫星遥感影像采用像元二分模型对伏牛山整体植被变化趋势进行研究，阿多等[14]对华北平原气候时空变化及其对植被覆盖度的影响进行了分析。裴志林等[15]对黄河上游植被覆盖度分布特征和影响因素进行了分析，发现气候类环境因素对植被覆盖度的影响较大。赵明伟等[16]对2001—2015年我国陆地植被覆盖度时空变化及驱动力进行了分析，发现对于不同地区影响其植被覆盖度变化的因素不尽相同。

伏牛山生态功能区对区域生态环境起着重要的调节作用，本文则在前人研究的基础上基于像元二分模型采用MODIS数据反演伏牛山地区的植被覆盖度，并采用一元回归、Hurst指数、马尔科夫矩阵等方法，分析研究该区域植被覆盖度的时空变化规律，探究气象因子对其的影响，为生态环境的监测提供依据并为下一步对研究区植被生态质量评估提供数据支撑。

1 研究区概况及数据来源

1.1 研究区概况

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

图1 研究区分布概况

1.2 数据来源及预处理

植被覆盖度的变化受气候影响较大。北半球植被覆盖度与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年月降水量和月均温变化

本研究采用的遥感数据来源于NASA的MODIS合成产品MODIS13A1，它是分辨率为500 m的16 d合成的L3级产品。时间覆盖范围是2000年4月至9月到2018年4月至9月，共19 a。河南省对应的正弦投影系统为h27v05，使用MRT软件进行投影转换，并将投影转换为Albers等面积割圆锥投影。基于ArcGIS工具根据研究区行政边界对数据做裁剪处理并进行NDVI的MVC月合成，得到每月对应的NDVI值。

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

2 研究方法

2.1 像元二分模型

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

NDVI=NDVIveg+NDVIsoil

(1)

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

NDVI=VFC×NDVIveg+(1-VFC)×NDVIsoil

(2)

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

(3)

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

表1 植被覆盖度等级划分

2.2 马尔科夫模型

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

(4)

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

2.3 一元线性回归分析法

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

(5)

式中：n是时间序列的长度(n=19)；VFCi是第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.010)，显著减少(0.010.05)。

2.4 Hurst指数

自然界中具有长程依赖性的时间序列是普遍存在的，并在水温，地球化学，气候，地质和地震等领域广泛运用[25]。而Hurst指数是描述长程依赖性的一种方法。Hurst指数可以反映出序列自相似性以及序列发展的相关强度[26]。基于重标极差(R/S)分析法计算得到的Hurst指数可以反映序列自相似性及序列发展的相关强度，并能够较好的判别时间序列数据的自相关性。若0

2.5 偏相关分析

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

(6)

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

3 结果与分析

3.1 植被覆盖度分布特点

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

3.2 植被覆盖度整体变化分析

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

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

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

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

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

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

3.3 植被覆盖度年际变化特征

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

3.4 植被覆盖度变化趋势

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

3.5 植被覆盖度与气象因子的相关性分析

由于气温和降水都会对研究区内的植被覆盖度的变化产生影响，将研究区内的植被覆盖度分别与月降水量和月均温进行偏相关分析，先选取代表年进分析，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年研究区森林地区植被覆盖度与降水量和温度的偏相关数

4 结论与讨论

(1) 研究区植被覆盖度整体较好。(2) 不同等级植被覆盖度类型转化的特点是由低水平覆盖度区向高水平覆盖度区的转化。(3) 自2000—2018年，研究区植被覆盖度变化趋势以增长为主，研究区内一半以上的地区其植被覆盖度呈现上升趋势。(4) 研究区的植被覆盖度变化的反向持续性要强于同向持续性，并且持续性明显，未来植被覆盖度的变化趋势预计以减少为主。(5) 研究区植被覆盖度的变化受气温影响较大。其中森林区域的植被覆盖度同样受气温的影响较大。这与阿多[14]等对华北平原气候时空变化及其对植被覆盖度的影响中的结论相一致即气温是影响森林生态区植被覆盖度的主导因子。本文与其均采用了较长的时间序列，可以较为充分地反映研究区植被覆盖度变化与气象因子之间的相关关系。张晓东[13]等在对伏牛山植被覆盖度的变化研究中，采用三期遥感影像，发现伏牛山地区的植被覆盖度变化和温度变化的关系不太明显，其推测是由于植被覆盖度的变化对气温的响应存有一定的滞后性导致的。由于采用的时间序列不同分析出的结果也不尽相同，综上对植被覆盖度与气象因子的分析采用长时间序列的遥感数据较为合理，可以降低气象因子的滞后性带来的影响。本研究没有考虑地形因素对研究区植被覆盖度的影响，研究地形因素对植被覆盖度的影响可作为以后的研究方向。