近10年汾河流域植被覆盖时空变化研究

韩丽君,李 晶

太原师范学院 地理科学学院,山西 晋中 030619

1 引言

1.1 选题背景与依据

土地利用/覆盖(LUCC)是全球变化研究的核心内容，也是区域生态评价的重要指标，在逐渐升温的环境变化和可持续发展研究中占据了重要地位。作为土地利用/覆盖最主要的表征，植被是最具有显示度的覆盖类型[1-11]。植被覆盖可作为生物地球化学过程、水文过程、环境变化的输入参数，在地学领域的研究中具有基础地位。

1.2 国内外研究现状

国内外植被覆盖提取的方法可分为两大类：一是地表实测，诸如目估法、样方法、仪器法等[12]。目估法是在实地或面对相片根据经验判断植被覆盖，此方法简单易行，但精度受人为因素影响较大。样方法是将研究区域划分为若干网格，收割网格内的全部植物，根据生物量和植被覆盖的比例系数推算植被覆盖，这种方法具有破坏性，且耗时耗力，多用于小区域点状研究对象的植被覆盖，难以大面积推广。仪器法通常使用相机拍摄植被冠层，利用图像处理软件提取植被覆盖，此方法测量结果具有较高精度，估算效率也较高，为目前最常用的实测方法。二是遥感监测法，诸如回归模型法、植被指数法、像元分解模型法等。回归模型法是利用回归分析得到空间外推公式，对影像波段计算出的植被指数或植被覆盖度进行空间或时间上的外延扩展，可得到更大范围或更多时期的植被覆盖，使用最广的为线性回归模型。如Larsson H.利用TM影像反演出阿拉伯地区的NDVI指数，并与植被覆盖度进行了线性回归，得到了精度较高的统计模型。Anatoly利用反演的NDVI、ARVI、GreenNDVI三种植被指数与小麦的植被覆盖进行了线性回归分析，得出了植被指数对于植被覆盖变化十分敏感的结论，并在一定程度上消除了大气的影响。从前人所做的研究来看，回归模型法存在一定的局限性，尽管对于局部区域的植被覆盖估算在满足一定效率的基础上具有较高的精度，但建立的回归模型只适用于特定的区域或特定的植被，不宜推广，同时所用数据多是Landsat、SPOT这种中高分辨率的影像，低分辨率的遥感影像由于混合信息较多不宜进行回归分析，同时，面对研究范围较大，时间序列较长的情况时会遇到数据获取的困难。像元分解模型模型法原理是将像元信息分解成多个组分(如土壤、水、植被等)，构建分解模型计算覆盖。植被指数法是建立植被指数与植被覆盖的关系，计算植被覆盖。

最重要的是，作为一种缓变自然现象，植被覆盖能够被遥感影像所记录。同时，由于环境条件是植被覆盖时空分布和变化的直接驱动力，植被对环境指标具有降低变化幅度等反馈的同时也能表征环境因子的变迁[13]。本研究通过植被覆盖的时空分布模式来定性与定量评价汾河流域生态环境，为区域的可持续发展提供科学依据。

2 模型与方法

2.1 估算模型

目前已经发展了许多利用遥感影像反演植被覆盖的模型，较为实用且简单的方法是利用植被指数估算植被覆盖。植被指数是遥感影像反演地表植被生长状况的一个指标[14]，常用的植被指数是归一化植被指数NDVI(Normalized Difference Vegetation Index)[15]，它利用了植被在红光波段的高吸收和近红外波段的高反射这个特征，计算公式为：

式中ρnir为近红外波段的反射率ρred为红光波段的反射率。NDVI值在[-1,1]之间[15]，正值表示有植被区，且植被状况越好值越高，负值或零表明像元为云、水、裸地等无植被区。

从植被指数到植被覆盖度的计算采用像元二分模型，它假设一个像元仅由两种地物(如植被区和无植被区)组成，影像记录的光谱信息是两种地物信息的线性加权和，各组分权重是其面积所占比例。

根据上述像元二分模型原理，像元信息P可分解为由植被贡献的信息Pveg和土壤、裸岩等无植被区贡献的信息Pbare[16]。

对于由植被区和无植被区构成的混合像元，植被区所占比例即是此像元植被覆盖度Cveg，与之对应的是非植被区覆盖：

式中，Pveg为完全有植被覆盖的纯像元所反映的光谱信息，Cveg*Pveg表示混合像元情况下植被所贡献的光谱信息；同理，Pbare为完全无植被覆盖的纯像元所反映的光谱信息，(1-Cveg)*Pbare表示混合像元情况下无植被区所贡献的光谱信息。

将上述公式进行方程组求解，可以得到：

像元NDVI值可视为有植被覆盖所贡献的NDVIveg和无植被覆盖所贡献的NDVIbare的组合，以此作为像元信息P代入公式，可以得到利用归一化植被指数计算植被覆盖度的公式：

式中，NDVIbare为完全无植被覆盖像元的NDVI值，NDVIveg为完全由植被覆盖像元的NDVI值。

利用这个模型计算植被覆盖度的关键是计算NDVIsoil和NDVIveg。这里有两种假设：

（1）当区域内可以近似取VFCmax=100%，VFCmin=0%。

公式（1）可变为：VFC=(NDVI-NDVImin)/(NDVImax-NDVImin)

NDVImax和NDVImin分别为区域内最大和最小的NDVI值。由于不可避免存在噪声，NDVImax和NDVImin一般取一定置信度范围内的最大值与最小值，置信度的取值主要根据影像实际情况来定。

（2）当区域内不能近似取VFCmax=100%，VFCmin=0%

当有实测数据的情况下，取实测数据中的植被覆盖度的最大值和最小值作为VFCmax和VFCmin，这两个实测数据对应影像的NDVI作为NDVImax和NDVImin。

当没有实测数据的情况下，取一定置信度范围内的NDVImax和NDVImin。VFCmax和VFCmin根据经验估算[17]。

一般情况下，可根据影像区域内NDVI分布，以0.5%的置信度的上下限阈值分别代表NDVIveg和NDVIbare。

基于像元二分模型的植被覆盖计算模型理论简单，所需参数容易获取，在实际应用中使用比较广泛。

2.2 分析方法

本研究在时间和空间两个维度对研究区植被覆盖发展变化及其与环境因子的关系进行分析：一是利用空间平均法分析植被覆盖时间变化及其与环境因子的关系，以发掘植被覆盖的时间分布模式；二是利用逐像元法分析植被覆盖空间变化及其与环境因子的关系，以发掘植被覆盖的空间分布模式。植被与影响因子的关系通过相关系数来研究。考虑到数据获取等因素的限制，本文计算了植被覆盖与气温、降水量的相关性。对于一个像元来说，10年际的16 d步长温度(或降水)会形成一个时间序列，同一时期的植被覆盖也会形成一个同等长度的时间序列，相关系数即可以得到每个像元植被覆盖与气象因子的相关性，公式为：。

式中i为16 d植被覆盖均值，j为16 d气温均值或降水量，相关系数为两个变量集的协方差除以它们的标准差乘积。

为定量化覆盖程度，等间距将VFC划分为5级反映覆盖程度，具体划分如表1。

表1 植被覆盖分级标准Table 1 Classification criteria for VFC

3 植被覆盖信息提取

3.1 研究区概况

汾河流域位于黄河中游，处于 34°N～38 °N、110°E～113°E 之间，行政区域包括忻州、太原、晋中、吕梁、临汾、运城6个地级市，40个县。命名河流汾河发源于山西忻州宁武县管涔山，南向跨越太原后在新绛县折向西，在河津市汇入黄河，全长709.9 km，流域面积39826 km2，是黄河的第二大支流。河源至太原为上游地区，多为山地、黄土丘陵，此段水土流失严重，是洪水、泥沙的来源地；太原至洪洞县为中游地区，进入晋中平原后坡降变缓，泥沙淤积严重。作为山西的母亲河，汾河在山西的自然、历史、文化地位举足轻重。汾河流域是山西人口密集、工农业集中的区域，人口占全省比例为39%，工业产值占全省46%，耕地面积占全省30%，农业产值占全省64%。历史时期的统计资料表明，中游地区的洪涝干旱频发。研究区如图1所示，右下角为研究区在山西的位置。

图1 研究区示意图Fig.1 Location of the studied areas

3.2 数据与处理

本研究利用VFC指数距平来估算植被覆盖，必须的数据来源于MODIS的NDVI产品。NDVI数据为美国航空航天局TERRA卫星搭载的中分辨率成像光谱仪(MODIS)的多时相合成产品。卫星每天都能获取同一区域的原始影像，为减少或避免气象原因导致的无数值区域，其产品数据采用一定的算法进行多天的合成。NDVI16 d合成产品采用了平均值算法，为展开研究工作，所用遥感影像来源于NASA LAADS网站(http://ladsweb.nascom.nasa.gov/data)，为2004～2013年MODIS产品数据MOD13A2，为1 km分辨率植被指数16 d合成产品。下载了图幅号为h26v5和h27v5的MODIS产品，数据为HDF格式，作为产品，已经进行了几何纠正和辐射纠正，通过镶嵌，投影，利用研究区矢量边界裁剪出研究所需的影像。均值合成后的NDVI数据在内部存在无数值区域，利用ArcGIS的Nibble功能修复坏点，形成内部连续的影像，统一为img格式。数据统一为六度分带的Transverse Mercator平面坐标投影，在这种区域尺度下，能最大程度减少区域内面积或角度属性的变形。所用参数如表2。

表2 Transverse Mercator投影参数Table 2 Projection parameters of Transverse Mercator

气象资料来源于国家气象中心(http://cdc.cma.gov.cn/home.do)，利用汾河流域及周边地区22个气象站2004～2013年降水、气温日值序列，以16 d累积量或平均值，使用反距离加权方法进行插值，分别得到230个时间序列的降水和气温面域化数据。数字高程模型DEM来源于ASTER GDEM(http://gdem.ersdac.jspacesystems.or.jp)。

4 植被覆盖分析

4.1 植被覆盖时空演化分析

利用2004～2013年16 d合成的MODISNDVI数据，根据上述估算模型，在 ArcGIS、Erdas、Envi、Matlab的协同工作下，结合流域边界矢量数据，实现汾河流域植被覆盖的遥感反演。图2显示了2004～2013年16 d步长的植被覆盖时空变化。

图2 2004～2013年16 d步长的植被覆盖时空变化Fig.2 Temporal and spatial variations of VFC for 16 d step length from 2004 to 2013

每年从年积日(Day of Year,DOY)第1 d到最后1 d依据16 d步长共有23个时相。综合10年际的时空分布，从DOY001到DOY081(对应1～3月)这一时期的植被覆盖变化不大，大部分区域的VFC值在0.4以下，只有少部分区域的植被覆盖比较可观，太原西部吕梁山区相对有一个聚集区VFC值在0.4以上。这一时期整体是两边的山区(吕梁山脉,太行山脉)VFC值较高，中间的河谷盆地值较低。到后期，VFC值在南部有一定的升高，临汾盆地相对明显。DOY097到DOY113这两个时相相对稳定，从DOY129开始，VFC值在流域边缘区域有明显的升高，太原盆地VFC值升高的幅度相对较低，在DOY129开始会形成一个明显的低值区域，并一直持续到DOY161。这一时期，流域两边的山区VFC值增加明显。在DOY177这一时相，流域北部，特别是太原盆地VFC值增加明显，这一时相北部的VFC值明显超过流域南部。从DOY193开始，流域VFC值开始大规模攀升，流域内值基本在0.6以上，两边山区的VFC值多在0.8以上，这一状态一直持续到DOY241和DOY257。这一时期过后，VFC值开始有较为明显的降低，在太原盆地和临汾盆地降速尤为显著。从DOY289开始，流域内VFC值降低到0.4以下，太原盆地的值最低。

4.2 植被覆盖与气象因子的关系

利用相关性来分析植被覆盖与气象因子的关系，两种气象因子对植被覆盖的影响力和控制力因地理位置的不同而异。图3显示了两种气象因子与每个像元16 d植被覆盖时间序列的相关性在空间上的分布。左边为降水与植被覆盖的空间相关性，降水与植被覆盖的相关性呈现出北部高、南部低的特点，整体上相关性偏低；右边为气温与植被覆盖的空间相关性，气温与植被覆盖的相关性呈现出四周高，中间低的分布，在介休县附近有一个相关性的低值区域，在太原北部偏西的区域相关性较高，普遍在0.8以上。

图3 植被覆盖与气象因子的关系Fig.3 The relationship between vegetation coverage and meteorological factors

图4 VFC、温度、降水的时间序列分布Fig.4 Time series distribution of VFC,temperature and precipitation

为发掘流域VFC的时间分布模式，将流域内VFC、气温、降水的均值进行联合对比分析。如图4所示，从上到下依次为VFC、气温、降水，红色曲线为最大值，绿色曲线为标准差，蓝色曲线为均值，蓝绿色曲线为最小值。坐标横轴为时间序列，每年23个时相，10年一共230个时间序列，横坐标为各指标值(为清晰显示做了偏置，观察变化趋势即可)。通过统计图与统计数字，发现VFC值与气温、降水在时间分布上有较好的一致性，但在时间上存在一个时相的步长差。VFC峰值出现在DOY209附近，气温和降水的峰值出现在DOY225附近。

4.3 植被覆盖与地形因子的关系

地形因子能影响植被的光热条件和水分状况，进而对植被覆盖产生作用。对区域内数字高程模型进行坡度分析和坡向分析，通过阈值分类，区分出坡度分类图和坡向分类图。

海拔分级标准按0～500 m、500～1000 m、1000～1500 m、1500～2000 m、＞-2000 m共5个级别。统计得到不同海拔的植被覆盖均值如图5所示。

图5 植被覆盖的海拔分布Fig.5 Statistics for vegetation cover altitudinal distribution

从图5可以看出，在海拔1500～2000 m这个区间的VFC值最高，在这个海拔之下，植被的VFC值随海拔的升高而升高，但幅度较小，植被分布随海拔高度的变化不大。在海拔大于等于2000 m这个区间的VFC值有明显的降低，说明这个海拔对对植被的影响较为显著，限制了植被的生长。通常来说，低海拔地区人类活动较为频繁，高海拔地区人类活动较少，低海拔区域的VFC值不高可能与此有关。

坡度分级标准为 0°～3°、3°～6°、6°～9°、9°～12°、＞=12°共5个类别，将坡度分类图与植被覆盖进行叠置分析，计算每个坡度范围内的植被覆盖均值，得到植被覆盖的坡度分布，如图6。

平均下来，0～3度范围内植被覆盖均值为0.3022，3～6度范围内植被覆盖均值为0.2199，6～9度范围内植被均值为0.1787，9～12度范围内植被覆盖均值为0.1533，大于等于12度地区植被覆盖均值为0.1431。坡度越大，VFC值显著降低，说明坡度是影响植被生长的一个主要因子。其原因在于坡度的大小与水分密切相关，坡度增加会导致水分流失严重，致使植被生长状况不佳。

对数字高程模型进行坡向分级，分级标准为135～225（阳坡）、90～135和225～270（半阳坡）、315～45（阴坡）、45～90和270～315（半阴坡）、无坡向共5类，如图7所示。将坡向图与植被覆盖进行叠置分析。计算每个坡向范围内的植被覆盖均值，得到植被覆盖的坡向分布，如图7。

图6 植被覆盖的坡度分布Fig.6 Statistics for vegetation cover slope distribution

图7 植被覆盖的坡向分布Fig.7 Statistics for vegetation cover slope distribution

从图7中可以看出，阴坡植被覆盖均值最高，半阴坡植被覆盖次之，半阳坡植被覆盖最低。结合上述坡度对于VFC值的影响，说明此区域降水是限制植被生长的一个重要因子，阴坡光照少，有利于水分的涵养，植被茂盛，半阳坡和阳坡因阳光照射导致水分蒸发较多，水分流失严重，植被反而不是很理想。

5 总结与展望

5.1 总结

本研究利用MODIS NDVI数据，基于像元二分模型对汾河流域2004～2013年10年际16 d步长的植被覆盖度进行了估算，对植被覆盖的时空演变过程进行了全程描述，并结合气象因子和地形因子给出了分析。研究得出了以下几点结论：

（1）植被覆盖在2004～2013年较为稳定，从DOY001到DOY081变化不大，基本为0.4以下，从DOY129开始，VFC值有明显的升高，在流域边缘的山区尤为显著，从DOY193开始，流域VFC值快速升高，流域内VFC值多在0.6以上，DOY257过后，VFC值回落到0.4以下。总体来说，中间的太原盆地和临汾盆地VFC值较低，两边的吕梁山区和太行山区VFC值较高；

（2）对植被覆盖与气象因子的空间分布关系与时间分布关系研究表明：空间上，植被覆盖与降水的相关性北部高、南部低，植被覆盖与气温的相关性边缘高、中间低；时间上，植被覆盖与降水、气温具有一致的波形，都是夏季高、冬季低，但在相位上存在差异，植被覆盖的峰值相对于降水与气温的峰值提前一个步长(16 d)；

（3）通过统计不同海拔、坡度、坡向的VFC值，可以得出流域植被集中分布在500 m以下，坡度小于3°的阴坡地区，符合实际情况。从结果来看，低海拔区域植被覆盖较为稳定，海拔大于等于2000 m后植被覆盖有一个明显的降低；

（4）植被覆盖是一种复杂地表过程，其成因受多方面因素影响，发展变化过程难以完全由自然因子来控制。

本文作为阶段性的成果，特别是在16 d这样一个高时间分辨率的步长下得到了流域植被覆盖，相对以往不存在时间分辨率或以季或年为时间分辨率的植被覆盖计算来说不失为一个有益的尝试，同时，本研究探索了植被覆盖变化的自然因子和驱动力，可为植被保持和恢复提供科学依据。

5.2 不足与展望

（1）本文统计了植被覆盖的海拔、坡度、坡向分布，但影响植被生长的自然因子远不仅此，地形对于水分的重分配而导致的植被生长不均影响显著。人为的因子和驱动力因数据获取困难而未进行相关研究，从本研究的初步结果来看，很有可能是比较重要的一环。可以将有关植被、人口密度、交通网络、水系分布等图层进行叠加分析得出人为影响，这一切期待更为开放的数据政策；

（2）本研究使用的VFC基于植被指数，MODIS影像获取的植被指数在高覆盖植被区的生长季节具备较高可靠性，在稀疏植被区(覆盖度小于15%)的非生长季节或无植被区可信度较低，这些区域或时期的植被覆盖估算需进一步的研究；

（3）植被覆盖对气象因子的反馈作用有待研究，如降水、温度适宜时，植被生物量在生长季会保持上升趋势，遥感反演的植被指数同步上升，小区域内的气温也应保持在一定范围，如遇干旱，植被受水分胁迫作用，一方面生长受限，植被指数降低，另一方面植被无足够水分供给蒸腾作用时，部分气孔关闭，致使小区域内温度升高；

（4）本文从宏观尺度研究了流域内的植被覆盖，但流域内人口空间分布极度不均，工业农业均有可观规模的分布，特别是农业及其附属的水利设施对植被覆盖的影响不可忽视，农作物的收割不可避免带来VFC值的突降，水利设施的灌溉条件又能极大影响植被覆盖，这两者都属于人类活动，涉及到人与自然的耦合关系。本文仅仅从自然属性出发，研究了植被覆盖与降水、温度的影响，忽视人类活动的影响会致使研究结果存在一定偏差，尤其是近10年城市化提速的情况下，建成区扩展对植被覆盖的影响应该是拉低作用，这在未来的研究中应予以考虑；

本研究是利用VFC对过去10年植被覆盖情况进行了研究，给出的只是相对的覆盖等级，但协同降水、气温、地形进行的分析，对既定事实的描述有助于对未来植被覆盖的恢复，在后续的研究中期待更长时间序列、更多因子的加入。

[1]钟莉娜,赵文武.基于NDVI的黄土高原植被覆盖变化特征分析中国水土保持科学,2013,11(5):57-62

[2]王兮之,陆 龙,梁钊雄.青海湟水流域植被覆盖度时空变化分析[J].水土保持研究,2010,17(6):172-176

[3]朴世龙,方精云.最近18年来中国植被覆盖的动态变化[J].第四纪研究,2001,21(4):294-302

[4]汪明霞,王卫东.植被覆盖度的提取方法研究综述[J].黄河水利职业技术学院学报,2013,25(2):23-26

[5]李 琳,谭炳香,冯秀兰.北京郊区植被覆盖变化动态遥感监测—以怀柔区为例[J].农业网络信息,2008(6):38-41

[6]姜 烨.孙建国,李 庆.基于像元二分模型的植被覆盖度遥感信息提取[J].科技信息,2013(10):168-170

[7]梁四海,陈 江,金晓媚,等.近21年青藏高原植被覆盖变化规律[J].地球科学进展,2007,22(1):33-40

[8]陈国荣,刘 健,施聪智,等.基于GIS和遥感技术的植被覆盖度反演模型[J].北华大学学报:自然科学版,2010,11(3):251-254

[9]秦 伟,朱清科,张学霞,等.植被覆盖度及其测算方法研究进展[J].西北农林科技大学学报:自然科学版,2006,34(9):163-170

[10]周兆叶,储少林,王志伟,等.基于NDVI的植被覆盖度的变化分析—以甘肃省张掖市甘州区为例[J].草叶科学,2008,25(12):23-29

[11]张学珍,朱金峰.1982～2006年中国东部植被覆盖度的变化[J].气候与环境研究,2013,18(3):365-374

[12]章文波,符素华,刘宝元.目估法测量植被覆盖度的精度分析[J].北京师范大学学报:自然科学版,2001,37(3):402-408

[13]龙晓闽,周忠发,张 会,等.基于NDVI像元二分模型植被覆盖度反演喀斯特石漠化研究——以贵州毕节鸭池示范区为例[J].安徽农业科学,2010(8):4184-4186

[14]李晓香,陈媛丽.基于NDVI的鹤壁市植被覆盖动态分析研究[J].河南科技2014(7):188

[15]詹起林.基于森林植被指数的快速生态评价[D].上海:上海师范大学,2012

[16]聂 勇,范建容,贺秀斌,等.水土流失遥感调查中植被信息提取与评价指标讨论[J].水土保持通报,2007,37(4):10-14

[17]蒋友严,韩 涛,王有恒,等.石羊河调水10a来民勤绿洲生态脆弱性变化[J].干旱区研究,2014,31(1):157-162