1996—2015年黄河源区植被覆盖度提取和时空变化分析

王俊奇， 王广军， 梁四海， 杜海波， 彭红明

（1.中国地质大学（北京）土地科学技术学院，北京 100083； 2.中国地质大学（北京）水资源与环境学院，北京 100083；3.内蒙古自治区煤田地质局勘测队，内蒙古 呼和浩特 010010； 4.青海省环境地质勘查局，青海 西宁 810007）

0 引言

植被是陆地生态系统的重要组成部分，其在保持水土、涵养水源等方面起着重要的作用。作为描述植被覆盖状况的重要参数，植被覆盖度的变化对于评价生态环境状况和土地沙化程度有着重要的指导意义［1-2］。在过去的30 多年，受全球气候变暖和人为活动的影响，素有地球“第三极”之称的青藏高原，特别是黄河源区出现草地不同程度的退化和土地沙化等一系列生态环境问题，已经成为全球气候变化研究的关键地区［3-5］。为保护黄河源生态环境，国家于2005年、2014年分别实施“三江源生态保护和建设一期、二期工程”［6-7］。2016年，正式启动三江源国家公园体制试点。特别是2019—2020年，习近平总书记四次实地考察黄河，并多次强调沿黄生态环境保护的重要性。从植被覆盖度变化角度开展黄河源区生态环境变化研究，对于深入理解青藏高原多年冻土区在气候变化和人类活动双重作用下的植被响应，以及为黄河源区生态环境保护和治理提供辅助决策具有重要的意义。

随着遥感技术的不断发展和进步，植被覆盖度的反演方法也逐渐的完善与成熟。常用的反演方法包括植被指数法、经验模型法和混合像元分解方法［8-10］。受限于遥感数据的空间分辨率，混合像元分解方法在包括黄河源区的三江源地区受到了广泛的关注和利用。混合像元分解方法包括线性光谱混合分解方法（Linear Spectral Mixture Analysis，LSMA）和非线性光谱混合分解方法（Nonlinear Spectral Mixture Analysis，NLSMA）［11-12］。当 前 以LSMA为基础的像元二分模型在三江源地区应用较多。如吴喜芳等［13］基于像元二分模型方法，对1982—2006 年NOAA/AVHRR NDVI 和2000—2012 年MODIS NDVI 数据进行植被覆盖度提取，结果发现黄河源区植被覆盖度呈现由东南向西北递减趋势，并且整体植被覆盖度呈增加趋势；张颖等［1］采用同样的方法对三江源植被覆盖度进行提取，结果发现植被覆盖度的整体趋势与上述基本一致；马天啸等［14］对黄河源2000—2010 年植被覆盖度进行时空变化分析，结果表明植被退化区平均坡度和海拔更低，说明人类活动可能对植被退化产生了影响；肖桐等［15］采用像元二分模型对2000—2011 年三江源植被覆盖度进行提取，结果表明黄河源地区草地覆盖度恢复特征明显。尽管像元二分模型提取的结果能够从整体上反映植被覆盖度的变化趋势，但像元二分模型其实质是两端元的LSMA 模型。由于在LSMA模型中，植被一般以一个固定端元表示，忽略了“同物异谱”现象，导致包括像元二分模型在内的LSMA 模型在面向多种植被覆盖区域适用性较差。而在LSMA基础上发展起来的MESMA分解模型由于考虑了“同物异谱”［16］，在面向多种植被覆盖或异质性区域适用性较好。如廖春华等［17］采用MESMA和LSMA模型分别与实测数据对比，结果发现MESMA 提取的结果与实测数据相关性更高。并且相较于LSMA 模型，MESMA 可以对植被进行分类，因此更适合在地表类型复杂多变的黄河源区进行土地信息的提取。安如等［18］采用MESMA 模型从可食牧草、毒杂草、裸地三种端元对2012 年三江源典型区域的草地退化进行研究，并与实测数据对比，精度高达84.2%；此外，安如等［19］、Li 等［20］采用优化后的MESMA 方法黄河源玛多县“黑土滩”进行遥感定量识别，相较于传统的目视解译有了较大的改善。

考虑到以上研究多采用黄河源地区2000—2012 年的数据，对于2012 年以来鲜有研究，且基于MESMA 理论基础的研究区主要集中在“黑土滩”等黄河源典型区域［21-22］，范围有限，也没有对毒杂草的面积数量变化及其对畜牧业的影响进行分析。因此本文拟采用MESMA 模型对1996 年、2000 年、2004 年、2006 年、2008 年、2013 年、2015 年共7 年的黄河源植被覆盖度进行提取，并采用转移矩阵和一元线性回归趋势分析方法对提取出的植被覆盖度进行后续分析，最后对1996 年和2015 年的毒杂草分布情况统计并分析其对畜牧业发展的影响，以期为黄河源区生态环境的研究提供新的视角。

1 研究区与数据

1.1 研究区概况

本研究区所指的黄河源区地理坐标为33°03′～35°29′N、95°52′～99°52′E，位于青藏高原东北部，研究区示意图见图1。高山、峡谷等地貌分布其中，平均海拔约4 200m，总面积约4.4 万km2，地势西北高东南低。黄河源区属于典型的大陆性高原气候，寒冷干燥，东南空气湿润，西北气候干燥，太阳辐射强，日温差大，降水多集中在6—9 月［5，23-24］。作为“世界屋脊”青藏高原的一部分，黄河源区的自身气候变化和生态系统均十分敏感。其主要植被类型包括小嵩草、藏嵩草、早熟禾和紫花针茅等可食牧草，以及棘豆、火绒草等毒杂草；土壤类型包括寒钙土，草毡土以及沼泽土等［25-26］。

图1 研究区示意图Fig.1 Location of the study area

1.2 数据来源与预处理

黄河源区受海拔和气候影响，植被生长季相对较短，非生长季植被枯萎，且存在霜冻和冰雪覆盖现象，因此本研究选取植被处于最大生长季的7—9月且受云覆盖影响较小的1996 年、2000 年、2004年、2006 年、2008 年、2013 年、2015 年7 个年度35 景Landsat 影像。由于部分年份的个别景影像受云影像较大，故采用相邻年份的影像替代（表1）。由于2013 年采用两景Landsat7-ETM+数据，因此首先要对其进行去条带处理。之后利用ENVI5.3 软件对数据进行辐射定标、FLAASH 大气校正、地形校正和拼接裁剪等预处理。地形校正中采用的DEM 数据为2015 年1 月6 号发布的ASTER GDEM V2。此外，在MESMA 模型建立感兴趣区时采用的辅助数据包括1∶100万植被图及1∶100万土壤图［27-28］。

表1 Landsat影像信息列表Table 1 Landsat image information list

2 研究方法

2.1 多端元混合像元分解模型

LSMA 模型的基本假设是，混合像元的光谱是由两种及两种以上的固定地物以线性组合的方式构成。其具体模型为

式中：rω为混合像元光谱反射率；fi为每一个端元的光谱反射率；ri为各端元的丰度值；εγ为残余误差；N为波段总数。

MESMA模型是对LSMA模型的一种改进，其原理为在每个像元中采用不同的光谱分解模型（即采用不同的端元组分）进行测试，选择模型误差最小的作为提取本像元的分解模型，从而避免了LSMA 的缺陷［17］，在实践应用中也更符合实际情况。MESMA模型中常用的端元包括植被、非植被、土壤、阴影四种端元，本研究中选择了植被-土壤-阴影三种端元组合。MESMA程序主要包括以下两个主要步骤［29-31］：①构建光谱库并选取最优端元，②端元模型的确定。

2.1.1 光谱库的建立和最优端元的选取

光谱库的构建通常包含参考端元和图像端元两部分，本研究在黄河源区1∶100 万植被分类图和土壤分类图的辅助支撑下，以图像端元作为光谱库的主要构建方式。本文中选取植被、土壤、阴影三种端元，其中植被包含小嵩草、藏嵩草、紫花针茅、早熟禾等可食牧草以及铺散亚菊、火绒草、棘豆等毒杂草。土壤包含寒钙土、草毡土、沼泽土等类型。由于在遥感影像成像过程中，地物阴影普遍存在，因此将阴影作为一种端元更符合实际情况［17］。最终通过对不同种类地物分别建立感兴趣区，共选取植被端元2 056 个，土壤端元258 个，并将这些端元构建为备选光谱库。

由于端元之间存在着相关性、可分性等一些不可避免的因素特点，需要通过计算端元的平均均方根误差（endmember average root mean square error ，EAR）和最小平均波谱角（minimum average spectral angle，MASA）［17，31］来选择最优端元。本文中利用ENVI插件“Viper Tools 2”工具（下载地址：www.vipertools.org）来计算EAR 和MASA。依据EAR 越小，端元越纯且MASA 越小，端元相关性越小的原则，且经实地确认，如表2所示，最终选取出小嵩草、嵩草等植被类型以及寒钙土、草毡土等土壤类型的端元光谱曲线。图2~3 所示为研究区主要的植被、土壤端元光谱曲线图。

2.1.2 端元模型的确定

本文选择植被、土壤、阴影三种端元，端元模型选择二端元和三端元模型进行测试。研究表明二端元模型适合自然地物景观的模拟分解，三端元模型适合人为扰动下的地物景观模拟分解［32］，综合考虑，选择三端元模型进行黄河源区MESMA分解。

根据理论研究，各组分在像元中的比例为0～1，由于误差存在，端元比例设置为1.05 为宜，因此设置组分值范围为-0.05～1.05，阴影最大比例为0.8，RMSE 和残差约束为0.025［17，29］。由2.1.1 节可以选取出最终的三种端元光谱库（包括阴影），而后通过三种光谱库中的光谱进行排列组合，并通过上述的约束条件选择最优的端元模型，将原始图像中的每个像元分解为植被、土壤、阴影，从而获得各组分的丰度图像。

表2 端元实地验证结果Table 2 Field verification of endmembers

图2 植被端元光谱曲线图Fig.2 Vegetation endmember spectral curve

图3 土壤端元光谱曲线图Fig.3 Soil endmember spectral curve

2.2 一元线性回归趋势分析

本文中采用一元线性回归趋势分析在空间上对植被变化进行分析，本方法在像元尺度上对植被覆盖度进行空间变化分析，根据每个像素植被覆盖度随时间的变化特征来反映整个空间的变化规律［33-34］，公式为

式中：n为监测年数；fvci为第i年植被覆盖度；slope为多年植被覆盖度线性拟合斜率。

拟合直线反映研究区多年内的植被覆盖度变化趋势，斜率为正，表示植被覆盖度增加，反之则表示减少；拟合斜率的绝对值越大，表示变化越显著。借鉴相关的研究成果［35］将变化趋势划为显著减少（slope≤-0.015）、轻度减少（-0.015＜slope≤-0.005）、基本稳定（-0.005＜slope≤0.005）、轻度增加（0.005＜slope≤0.015）、显著增加（slope＞0.015）五个等级。

3 结果与分析

3.1 植被覆盖度空间分布特征

根据植被覆盖度提取结果并结合研究区的植被盖度实际情况，本文将植被覆盖度划分为5 个等级。即：植被覆盖度在0～10%为裸地，植被覆盖度在10%～30%为低覆盖度区域，植被覆盖度在30%～45%为中低覆盖度区域，植被覆盖度在45%～60%为中覆盖度区域，植被覆盖度在60%～100%为高覆盖度区域。运用MESMA模型对黄河源1996年、2000年、2004年、2006年、2008年、2013年、2015年共7年的数据提取植被覆盖度，并通过ArcGIS对植被覆盖度进行分级统计。由于居民区及荒漠化严重的北部等地区植被稀少，在本文的方法中只存在土壤一个端元，因此此处的植被覆盖度按裸地统计，详见图4（a～g）。研究区整体植被覆盖度东南高西北低，北部地区荒漠化严重。1996—2015年的植被覆盖度统计结果见表3。各年份的植被覆盖度统计说明具体如下：

图4 1996—2015年黄河源区植被覆盖度时间序列分布Fig.4 Time series of vegetation coverage in the source region of the Yellow River from 1996 to 2015

表3 黄河源1996—2015年植被覆盖度面积及比例Table 3 Area and proportion of vegetation coverage in the source region of the Yellow River from 1996 to 2015

（1）1996—2015 年，研究区裸地在2013 年占比最低，为10.60%，2004 年占比最高，为24.93%；其他年份裸地占比与2013 年略有浮动，但差异不大。研究区裸地主要分布在北部干旱地区及一些高海拔地区。

（2）研究区低覆盖、中低覆盖及中覆盖的区域面积占比基本保持一致，其中低覆盖面积占比在17%～20%，中低覆盖和中覆盖面积占比在12%～14%，且以上三种植被覆盖度区域主要集中在西部地区、扎陵湖和鄂陵湖的南部、鄂陵湖的东部及中部部分区域。

（3）研究区高覆盖区域面积占比结果显示，2004年面积占比最低，为28.65%，1996年、2000年、2006 年及2008 年面积占比约为40%，而2013 年、2015 年两年的面积占比相对较高，均在45%左右浮动，这也显示出近年来黄河源区植被覆盖度有所改善。研究区高植被覆盖区域主要分布在中南部、东部及东南部区域。

3.2 不同植被覆盖度动态变化分析

本文利用ArcGIS 软件对1996—2004 年、2004—2015 年两个时间段的植被覆盖度分别建立转移矩阵，获取植被覆盖度的动态变化特征，并对1996年、2004 年、2015 年3 年的不同土地覆盖类型制作面积比例图，如图5所示。转移矩阵详见表4～5。由上述两表并结合图5 分析可得，1996—2004 年裸地整体增加10.24%，未转化面积占比高达77.22%。低覆盖类型、中低覆盖类型和中覆盖类型整体基本保持不变，且转化主要向更低的植被覆盖度类型，分别占比为46.05%、61.38%、61.89%。而高覆盖类型整体降低12.75%，且未转化面积仅占61.05%。因此1996—2004 年植被覆盖度总体呈退化趋势。而2004—2015 年，裸地整体降低12.91%，主要转化为低覆盖类型，占比为39.09%。低覆盖类型、中低覆盖类型及中覆盖类型基本保持不变，且均主要向更高植被覆盖度类型转变，分别占比为62.87%、67.91%、72.91%。 而高覆盖类型整体增加16.18%，且未转化面积占比高达92.23%。综合以上分析，说明2004—2015年黄河源区整体植被盖度呈改善趋势。这也同时表明国家于2005年、2014年实施的“三江源生态保护和建设一期、二期工程”收到成效。

图5 1996年、2004年、2015年黄河源区不同土地覆盖类型的面积比例Fig.5 Area ratio of different land cover types in the source region of the Yellow River of 1996、2004、2015

3.3 植被覆盖度变化趋势分析

本文通过一元线性回归趋势分析方法（slope）对1996—2015年的植被覆盖度进行变化趋势分析，见图6，并对五种趋势进行统计，见表6。

由表6 可知，植被覆盖度呈增加趋势的区域面积占57.25%，其中显著增加占40.08%，轻度增加区域面积占17.17%，显著增加区域主要分布在中部、中东部、环鄂陵湖区域及北部部分区域，轻度增加区域则相对散乱，分布在研究区的各个部分。考虑到北部区域以裸地及低覆盖植被为主，而北部植被覆盖度呈增加趋势，说明北部地区生态环境改善较为明显；植被覆盖度呈下降趋势的区域占26.73%，显著降低区域占15.47%，主要分布在黄河源区的西部、扎陵湖北部、源区东北部及东南部部分区域，分析其原因，主要是因为黄河源头及一些河谷地带、环湖地区受人类影响较大，且东南部地区海拔较低受过度放牧影响较大以及扎陵湖北部的牧场消失等。轻微减少的区域占11.26%，分布相对散乱。基本不变的区域主要集中在北部及东部沙漠地区，占16.02%。因此黄河源区的治理范围应主要集中在常年不长植被地区及植被覆盖度下降区域，这也是生态监测的重点地区。

表4 黄河源1996—2004年植被覆盖度转移矩阵Table 4 Area transfer matrix of vegetation coverage in the source region of the Yellow River from 1996 to 2004

表5 黄河源2004—2015年植被覆盖度转移矩阵Table 5 Area transfer matrix of vegetation coverage in the source region of the Yellow River from 2004 to 2015

图6 黄河源1996—2015年植被覆盖度变化趋势Fig.6 Vegetation coverage tendency in the source region of the Yellow River from 1996 to 2015

表6 黄河源1996—2015年植被覆盖度变化趋势面积及比例Table 6 Area and proportion of vegetation coverage change in the source region of the Yellow River from 1996 to 2015

3.4 毒杂草分析

3.4.1 毒杂草分布

图7 1996年和2015年火绒草、棘豆、铺散亚菊分类Fig.7 Classification of Edelweiss，Oxytropis，and Ajania khartensis in 1996 and 2015

图8 1996年和2015年毒杂草植被覆盖度Fig.8 Poisonous weeds vegetation coverage in 1996 and 2015

本文采用三端元模型对源区铺散亚菊、火绒草、棘豆等主要的毒杂草类型进行提取，分析了1996年、2015年两年的三种毒杂草的分布概况并对毒杂草的总面积进行统计，详见图7（a、b）、图8（a、b）。从图7分析可得，黄河源区的毒杂草主要以火绒草、棘豆为主，火绒草主要分布在研究区南侧，棘豆则主要分布在北侧，铺散亚菊分布地区则相对分散。这与毒杂草的生长环境息息相关。火绒草较宜疏松土壤环境，棘豆则生长于沙地、山坡草地，而相较于整个研究区，北部主要以沙地为主，土壤较为贫瘠，南部地区土壤则相对疏松，土壤有机质含量较高。对图8 的毒杂草面积统计并结合表6 分析可得，1996—2015 年黄河源区总体植被覆盖度呈增加趋势，说明生态环境有所改善，但毒杂草的面积由1996 年的16 060 km2到2015 年的22 942 km2，整体增加6 882 km2，说明部分植被覆盖度增加的区域伴随着毒杂草的增加，因此对于部分黄河源区的畜牧业来讲是不利因素。

3.4.2 毒杂草提取结果验证

本文重点对MESMA提取出的植被类型从以下两方面进行了验证：

（1）结合实地考察数据，对18 个地点的植被类型提取结果进行了验证，验证结果如表7 所示。结果表明，MESMA 反演的18 个地点的植被类型数据，其中与实地考察结果一致的有13个地点，另有5个地点的提取结果与实地不一致。植被类型提取结果的正确率为72.22%。

（2）考虑到实地考察中有些地方难以到达，本文同时还利用2013 年8 月16 日获取的HJ-1A 高光谱遥感数据，基于堆栈式稀疏自编码器方法［36］，对37.76 km2（34°58′58″～35°02′12″N，97°52′13″～97°52′22″E）的无云覆盖区域进行了植被分类试验，将高光谱遥感数据的分类结果［图9（a）］与本文MESMA反演的植被类型提取结果［图9（b）］进行了对比。本文提取的毒杂草分布状况除个别区域外，总体空间分布状况与高光谱提取结果近似一致。在验证区内，MESMA 模型提取的毒杂草面积为18.56 km2，占比49.2%；堆栈式稀疏自编码器方法提取的毒杂草面积为15.22 km2，占比为40.3%；而两者提取出的公共毒杂草面积为9.35 km2。如果以HJ-1A高光谱数据分类结果为准的话，本文MESMA模型提取出的毒杂草的正确率为61.43%。但考虑到HJ-1A 高光谱遥感数据采用的是硬分类方法，对于毒杂草覆盖度较低的区域无法进行有效分类，以及HJ-1A 高光谱遥感数据空间分辨率较低的事实（HJ-1A HSI数据空间分辨率为100 m），认为本文的提取结果精度应高于61.43%，研究结论是可信的。

表7 端元一致性确认结果Table 7 Endmember consistency validation results

图9 毒杂草对比验证结果Fig.9 Comparison results of poisonous weeds：results of hyperspectral classification（a）and results of MESMA（b）

4 结论与展望

本文利用黄河源区七期陆地卫星影像，辅以1∶100 万植被图和土壤图，基于MESMA 模型完成了黄河源区植被覆盖度的提取，结论主要如下：

（1）研究区以中低植被覆盖度为主，且植被覆盖度由东南向西北呈递减趋势，这与研究区的海拔和气候条件相一致。

（2）研究区在1996—2015 年植被覆盖度总体呈改善趋势。其中1996—2004 年生态环境呈退化趋势，2004—2015年研究区生态环境有所改善。

（3）虽然研究区在1996—2015 年植被覆盖度总体呈改善趋势，但毒杂草的覆盖面积也在随之增加，不利于当地畜牧业的发展。

MESMA 模型对不同的像元采用不同的端元模型，提高了植被覆盖度的反演精度，且MESMA模型同时可以逐像素获取像元内的植被类型，从而提取出研究区的毒杂草。但由于条件限制，本文仅采用图像光谱，如何将实测光谱与图像光谱（包括星载高光谱影像）相结合，以提高植被类型提取的精度，也是未来考虑的方向。

本文所采用的MESMA模型提取出的黄河源区植被覆盖度变化在空间分布上与吴喜芳等［13］、张颖等［1］的研究基本保持一致，总体呈现出由东南向西北逐级递减的趋势。在时间分布上，1996—2004 年植被覆盖度变化趋势与康悦等［4］的研究一致，呈退化趋势；而在1996—2015 年与吴喜芳等［13］在相同时段内的研究一致，总体呈增加趋势。因为已有的研究已经证明了MESMA进行植被覆盖度提取的有效性［37-38］，所以本文主要是侧重进行毒杂草类型的验证，主要通过与实地考察数据和HJ-1A HSI 高光谱影像提取结果对比进行验证说明，且认为本文的植被类型提取结果精度应高于61.43%，研究结论是可信的。本文的研究结果反映了毒杂草在黄河源区广泛分布且逐步扩大的事实，畜牧业管理部门应对此加以重视。