基于地形梯度的羌塘国家级自然保护区土地利用时空演变特征及驱动因素

Spatiotemporal Evolution Characteristics and Driving Factors of Land Use in Chang Tang National Nature Reserve Based on Terrain Gradient

摘要

羌塘国家级自然保护区是青藏高原野生动物最集中、最具高原生态特征的生态地理单元，分析羌塘保护区的土地利用变化具有重要的现实意义。文章利用长时间序列高精度土地利用数据，运用地形梯度分级、PLUS模型和地理探测器，系统分析了保护区土地利用的时空演变特征、空间分异和转型驱动因素。结果表明：（1）1990—2021年研究区土地利用以草地和未利用地为主，二者约占95.13%；其次为冰川积雪和水域，二者约占4.87%。整体呈现水域扩张，其余地类缩减的趋势。（2）草地集中分布在海拔4 939~5 405m、坡度4.17°~25.93°的区域；水域高度集中分布在海拔4 939 m、坡度11.923°的平坦地区；冰川积雪主要分布在海拔5 597 m、坡度16.09°的阴坡；未利用地集中分布在海拔4 818 m、坡度8.05°的阳坡。（3）NDVI为土地利用空间分异的主导因子（q值=0.40），其次为年均温、海拔、距湖泊距离。（4）距居民点距离和海拔分别为草地和冰川积雪面积变化的主导因子，贡献度为0.14和0.40；距湖泊距离为水域和未利用地面积变化的主导因子，贡献度分别为0.26和0.14。综上，羌塘自然保护区需贯彻科学的国土空间规划，促进区域生态环境的可持续发展。

Abstract

The Chang Tang National Nature Reserve​ on the Tibetan Plateau​ is one of the most representative ecogeographic units with the highest concentrations of wildlife and the most distinctive plateau ecological characteristics​ in China. Consequently, examining land-use change​ in the Chang Tang National Nature Reserve had clear and important practical significance. The spatiotemporal dynamics​, spatial heterogeneity of land use and driving factors​ of land-use transitions in the Chang Tang National Nature Reserve were systematically analyzed based on the long time series of high-resolution land-use data​ by using topographic gradient stratification, PLUS model and geographical detectors. The findings were as follows: (1)​ Grassland and unused land dominated the land use types in the study area from year 1990 to 2021, jointly accounting for approximately 95.13%, followed by glacier–snow​ cover, and water body area, accounting for 4.87%. A clear trend of water-body expansion​ and contraction of other land types​ emerged. (2)​ Grassland was concentrated at elevations of 4 939\sim5 405 m​ with slope degree of 4.17°\sim25.93°; water body were highly concentrated at elevations 4 939 m​ with slope degree 11.923°​on flat terrain; glacier–snow cover was primarily distributed above 5 597 m​ with slope degree 16.09°​on northerly aspect; unused land was concentrated above 4 818 m​ with slope degree 8.05°​ on southerly aspect. (3)​NDVI​ was the leading factor shaping the spatial heterogeneity of land use (q=0.40), followed by mean annual temperature, elevation, and distance to lakes. (4)​ Distance to settlements​ and elevation​ were the primary drivers for changes in areas of grassland​ and glacier–snow​cover with contribution degree of 0.14​ and 0.40, respectively; distance to lakes​ was the primary driver of changes in area of water body​ and unused land with contribution degree of 0.26​ and 0.14, respectively.

引 言

羌塘国家级自然保护区位于青藏高原腹地，是全球海拔最高、面积最大的自然保护区之一，覆盖西藏自治区那曲市与阿里地区的核心生态区域^[1]。作为青藏高原高寒生态系统的典型代表^[2]，该区域不仅是藏羚羊、野牦牛等珍稀濒危野生动物的集中栖息地^[3]，更因其独特的高原气候与地形条件形成了以草地、冰川积雪、未利用地为主体的特殊土地利用格局^[4]。然而，受全球气候变化与人类活动干扰的双重影响^[5]，羌塘保护区的土地利用正经历显著变化：既有冰川退缩、湖泊扩张等自然过程驱动的地类转型^[6]，也存在因放牧压力、基础设施建设等引发的土地利用扰动^[7]。在此背景下，系统分析该区域土地利用的时空演变特征及其驱动机制，不仅能为保护区国土空间规划的精准制定提供科学依据，更有助于理解高寒生态脆弱区人地关系演化的规律，对维护青藏高原生态稳定性具有重要的理论与实践意义。

土地利用/覆被变化（LUCC）是全球环境变化研究的核心内容之一^[8]。20世纪90年代以来，国际地圈-生物圈计划（IGBP）与全球环境变化人文因素计划（IHDP）联合推动的LUCC研究，已从早期的全球或大尺度宏观分析，逐步转向区域尺度的精细化探究^[9]。针对高寒生态区，学者们围绕青藏高原开展了大量工作：一方面，聚焦土地利用类型的数量与空间转型特征，指出草地长期是高寒生态区的主导地类，但受气候变化与超载放牧影响，部分地区存在草地退化与沙化趋势^[10-12]；另一方面，关注驱动因素的作用机制，指出自然因子（如海拔、温度、降水）通过限制植被生长与地表过程塑造了土地利用的基本格局^[11]，而人口密度、畜牧业规模等社会经济因子则通过直接或间接方式影响地类转换强度^[10]。值得注意的是，地形作为连接自然环境与人类活动的关键介质，其对土地利用分布与变化的调控作用逐渐受到重视。已有研究表明，海拔通过影响水热组合条件形成垂直地带性分异，进而约束不同地类的适宜分布范围^[13]；坡度与坡向则通过改变光照、水分再分配及人类可进入性，进一步细化土地利用的空间格局^[14]。然而，现有针对羌塘保护区的相关研究多集中于单一地类的变化分析^[15]，缺乏从地形梯度视角系统分析多类型土地利用（草地、水域、冰川积雪、未利用地等）的时空演变规律及其驱动机制的研究，尤其对地形、自然因子、人类活动多维驱动力的交互作用关注不足。

本文以羌塘保护区为研究对象，以1990—2021年长时间序列高精度土地利用数据为基础，采用地形梯度分级方法，将研究区按照海拔、坡度、坡向划分为多级梯度单元，结合PLUS模型模拟地类转换过程，并利用地理探测器量化自然与经济驱动因子的贡献度及交互效应。该方法能够突破传统研究中对地形因子静态背景的简单控制，能够有效揭示不同地形条件下土地利用变化的异质性特征，为理解高寒生态区地形、土地利用、驱动因子的耦合机制提供新的分析视角。

1 研究区概况

羌塘自然保护区位于西藏自治区西北部（图1），涵盖那曲市西北部和阿里地区东北部^[16]，是仅次于格陵兰国家公园的全球第二大陆地自然保护区^[1]，也是地球上最后一个基本未受干扰的牧场生态系统^[17]。该区地貌以低山、丘陵与湖泊为主，是著名的高海拔湖泊集中区，平均海拔在5 000 m以上，被称为世界屋脊的屋脊，是全球平均海拔最高的自然保护区。羌塘保护区空气稀薄，气候寒冷而干燥，年均温在0℃以下，年均降水约154.9 mm，冬春多大风^[1]。羌塘自然保护区设立于1993年，2000年晋升为国家级自然保护区，以藏羚羊、藏野驴、野牦牛和藏棕熊等濒危野生动物为保护对象^[17]。该地区人烟稀少，大部分区域直至20世纪后半叶才开始有人居住^[18]。

2 数据来源与研究方法

2.1 数据来源

①逐年土地利用数据（1900—2021年）来源于中国土地覆盖年度数据集（http:// https://zenodo.org/record/5210928），分辨率为30 m，总体精度为79.31%^[19]，满足研究需求。②年均温、年均降水、人口密度、GDP、NDVI数据来源于中国科学院资源环境科学数据中心（http://www.resdc.cn/）。③DEM数据来源于地理空间数据云（http://www.gscloud.cn/），分辨率为30 m，并由此处理得到坡度和坡向。④公路、居民点、水域、水系数据来源于国家基础地理信息中心（https://www.webmap.cn/）。

2.2 地形数据处理及分布指数

选取海拔、坡度和坡向3个地形要素分析羌塘保护区土地利用的地形梯度变化。其中，海拔和坡度通过自然断点法^[20]进行分级，坡向依据前人研究成果^[21]进行分级（表1）。

为消除各地形特征面积大小差异在土地利用分布中的影响，采用标准化、无量纲的分布指数^[22]来描述土地利用在不同地形梯度上的分布情况，公式如下：

P=\left(S_{ie}/S_i\right)/\left(S_e/S\right)（1）

式中，P为分布指数；S为研究区总面积；S_e为研究区e地形要素的土地利用总面积；S_i为i类土地利用总面积；S_{ie}为e地形要素下i类土地利用面积。P1时，即该地类在此地形区间上为优势位^[23]。

2.3 地理探测器

地理探测器可探索要素的空间分异特征^[24]，采用因子探测器和交互探测器对羌塘自然保护区土地利用空间分异的驱动因子及其相互作用进行探测，公式如下：

q=\mathrm{1}-\frac{\mathrm{1}}{N{\mathrm{\sigma }}^{\mathrm{2}}}{\displaystyle \sum _{h=\mathrm{1}}^L}{N}_h{{\mathrm{\sigma }}_h}^{\mathrm{2}}（2）

式中，q为某因子对土地利用空间分异的影响力；h为变量Y或因子X的分层/分区；{{\mathrm{\sigma }}_h}^{\mathrm{2}}和{\mathrm{\sigma }}^{\mathrm{2}}分别为评价单位和全区Y值的方差；N_h和N分别为h区和全区单元数。

2.4 土地利用扩张驱动模型

PLUS模型中的用地扩张分析策略采用随机森林分类（RFC）算法逐一对各地类扩张和驱动因素进行挖掘，可获取特定时段各地类的发展概率及驱动因素对各地类扩张的贡献^[25]。RFC算法能处理高维数据及变量之间的多重共线性，确定k类土地利用在单元格i上出现的概率P_{i,k}^d，其演算公式^{[25, 26]}如下：

P_{i,k}^d\left(x\right)=\frac{{\sum }_{n=\mathrm{1}}^{M}I\left(h_n\left(x\right)=d\right)}{M}（3）

式中，d取0或1，取1表明由其他地类转为k类，0则表示土地利用类型转变成了除k外的其他地类；x为多个驱动因素组成的向量；I为决策树集的指示性函数；h_n\left(x\right)为向量x的第n个决策树的预测类型；M为决策树总数。

3 结果与分析

3.1 土地利用的时空演变特征

羌塘自然保护区土地利用具有显著的时空分异性（图2和图3）。在空间上，羌塘自然保护区内草地在中部唐古拉山脉南麓地势低缓地带呈连片状、北麓呈狭窄带状分布（图2）；水域多以湖泊形式集中于保护区东北部、东南部和西部；冰川积雪高度集中于中部唐古拉山脉、西北部昆仑山脉和东北部可可西里山脉的高海拔山峰区域；未利用地以广阔带状集中于北部昆仑山脉南麓。在时间上，1990—2021年保护区内土地利用类型以草地为主导（70.27%），未利用地占比也较大（24.86%），其次为水域（3.32%）和冰川积雪（1.55%），耕地、灌木和建设用地极少（图3）。从面积变化来看，1990—2010年草地呈减少趋势，2010—2021年则呈增长趋势，未利用地刚好相反，水域持续增长，冰川积雪持续减少；1990—2021年草地、未利用地和冰川积雪面积分别退化1.46%（2 859.98 km²）、2.00%（1 285.73 km²）和27.16%（1 411.36 km²），水域扩张87.30%（5 558.28 km²），其余地类面积变化极少，表明1990—2021年羌塘保护区土地利用总体呈现草地、水域和未利用地面积缩减，水域则扩张的变化规律。

从土地利用转移来看，1990—2021年羌塘保护区主要土地利用类型间转换频繁（图4），1990—2005年未利用地和水域在大量草地和冰川积雪的转入下，面积增长迅速，其中有14 835.94 km²草地和721.81 km²冰川积雪转为未利用地，集中于保护区北部和西部高寒区（图3）；此外，草地和冰川积雪转入水域的面积分别为1 409.98 km²和195.67 km²，集中于保护区东南区湖泊连片带。2005—2021年草地和水域在未利用地和冰川积雪转入下，面积迎来扩张，分别有12 830.49 km²未利用地和46.22 km²冰川积雪转为草地，集中于保护区西部广袤的山麓地带；水域扩张面积集中在保护区北部和西部湖泊集中带，分别有1 723.09 km²未利用地和107.78 km²冰川积雪转入。总体而言，羌塘保护区土地利用在全球气候变暖等自然因素、过度放牧和退牧还草等自经济因素的双重影响下，冰川积雪持续萎缩，水域面积不断扩张，未利用地和草地在2010年前后受政策影响，面积增长趋势截然不同。

3.2 土地利用的地形梯度分布和变化

羌塘自然保护区地势大致呈西北高东南低的态势（图1），受东西走向的昆仑山、唐古拉山和可可西里山等高大山脉影响，北部和南部地形较为平坦。85.54%的区域海拔介于4 818～5 405 m，约92%的区域坡度25.93°，海拔和坡度较低区域为保护区北部、东南部和西部湖泊集中地区。除以湖泊为主的平面（1.87%）较少外，各坡向面积占12%左右（表1）。

在地形梯度分布方面，以2021年为代表年，得到羌塘保护区主要土地利用类型的地形梯度分布情况（图5和图6）。水域面积随海拔上升呈先增后减趋势，随坡度和坡向则呈递减趋势。其分布集中于第Ⅱ级海拔、第Ⅰ级坡度和坡向内，分别分布为42.64%（5 084.52 km²）、58.09%（6 926.36 km²）和41.87 %（4 992.42 km²），第Ⅰ\simⅡ级海拔、第Ⅰ级坡度和坡向为其分布优势带，其中第Ⅰ级坡向分布指数高达22.33。冰川积雪面积随海拔上升大致递增，第Ⅸ级内分布有52.92%（2 003.00 km²），分布指数高达54.97。第Ⅶ\simⅨ级海拔为其分布优势位，集中了98.10%的冰川积雪；冰川积雪面积在坡度和坡向尺度呈先增后降的变化趋势，第Ⅴ级坡度分布14.06%（530.99 km²），除平面仅分布0.02%（0.85km²）外，各级坡向分布面积占10.41%（393.05 km²）到14.77%（557.86 km²），第Ⅴ\simⅨ级坡度和第Ⅱ\simⅤ、Ⅸ级坡向为其分布优势位，阴坡内分布55.23%（2 086.24 m²）。未利用地面积随海拔、坡度和坡向增大总体呈先增后降的变化趋势，其中第Ⅲ级海拔、第Ⅱ级坡度和第Ⅶ级坡向内分布最多，分别占23.60%（14 894.86 km²）、22.24%（14 035.01 km²）和13.24%（8 357.05 km²）。第Ⅱ~Ⅲ级、第Ⅶ~Ⅷ级海拔，第Ⅱ和Ⅵ\simⅨ级坡度，第Ⅱ和Ⅴ\simⅧ级坡向为其分布优势位。草地面积随海拔和坡度上升呈先增长后减少的趋势，其中第Ⅲ级海拔和第Ⅱ级坡度内分布最多，面积分别占22.69%（43 900.44 km²）和21.64%（41 864.57 km²），第Ⅲ\simⅥ级海拔、第Ⅱ\simⅥ级坡度为其分布的优势带；坡向尺度上除在平面仅分布0.04%（78.22 km²）外，其余各级内草地占比均在12%（23 200 km²）左右；第Ⅱ\simⅨ级坡向为其优势带。

在地形梯度分布变化方面，选取1990年、2005年和2021年为代表年，进一步分析羌塘保护区主要土地利用类型在地形梯度的时间变化（表2）。1990—2005年，除Ⅸ级海拔内草地面积增长145.22%（192.53 km²）外，其余各级海拔、坡度和坡向上的草地面积均不同程度减少，其中Ⅷ级海拔和坡度、Ⅰ级坡向内减少尤为显著，分别减少29.51%（654.14 km²）、10.90%（477.42 km²）和82.10%（443.39 km²）。2005—2021年，仅在Ⅰ\simⅡ级、Ⅸ级海拔和Ⅰ级坡向内草地面积有所缩减，降幅分别为1.78%（143.94 km²）、0.05%（14.35 km²）、3.34%（10.85 km²）和19.10%（18.47 km²），其余各尺度内草地呈增长趋势，其中Ⅷ级海拔、Ⅸ级坡度和坡向分别增长19.19%（299.91 km²）、11.85%（144.01 km²）和2.89%（678.75 km²）。未利用地面积变化情况与草地相反。1990—2005年，未利用土地在Ⅰ级、Ⅸ级海拔和Ⅰ级坡向内分别减少了2.11%（54.23 km²）、3.60%（10.60 km²）和84.81%（318.49 km²），其余各尺度均呈增长趋势，Ⅶ级海拔、Ⅸ级坡度和Ⅲ级坡向增长显著，分别增长27.70%（1 641.79 km²）、37.35%（259.20 km²）和9.56%（740.57 km²）。2005—2021年，仅Ⅸ级海拔的未利用地增加3.22%（9.14 km²），其余各尺度内均不同程度缩减，其中Ⅰ级海拔和坡向减少22.22%（560.40 km²）和43.50%（24.81 km²）、Ⅷ级坡度减少13.02%（315.70 km²）。水域和冰川积雪面积变化趋势也恰好相反，其中，水域除Ⅶ\simⅧ级海拔在1990—2005年减少11.71%（1.98 km²）和6.93%（1.56 km²）、2005—2021年减少10.75%（km²）和23.80%（4.97 km²）外，其余各尺度不断增长。其中，1990—2005年Ⅵ级海拔和坡向、Ⅴ级坡度水域增长191.47%（25.73 km²）、103.29%（263.46 km²）和136.24%（126.51 km²）；2005—2021年Ⅱ级海拔、Ⅳ级坡度和Ⅵ级坡向分别增长33.75%（1 283.12 km²）、75.00%（164.52 km²）和59.94%（310.77 km²）。2005—2021年，冰川积雪仅在Ⅸ级海拔略微增长0.01%（0.26 km²）外，其余各尺度内均呈退缩趋势。其中，1990—2005年Ⅰ级海拔、坡度和坡向分别减少88.69%（69.78 km²）、34.80%（142.16 km²）和79.90%（70.12 km²）；2005—2021年Ⅰ级海拔、坡度和坡向也减少显著，分别减少82.94%（7.38 km²）、18.45%（49.14 km²）和95.19%（16.79 km²）。

3.3 土地利用的空间分异驱动因素

空间分异性为自然和经济因子复杂耦合作用过程的空间表现，选取年降水（X1）、年均温（X2）、人口密度（X3）、NDVI（X4）、GDP（X5）、距湖泊距离（X6）、距河流距离（X7）、距居民点距离（X8）、距公路距离（X9）、海拔（X10）、坡度（X11）和坡向（X12）共12个自然和经济因子进行探测，羌塘保护区土地利用空间分异的驱动因子解释力由高到低依次为：X4X8X1X2X10X9X3X5X6X11X7X12（图7）。其中，X4的解释力最为显著，达0.41；其次为X1、X2和X8，解释力在0.171以上，是影响保护区土地利用空间分异的重要因素；X3、X9和X10的解释力在0.04以上，是影响保护区土地利用空间分异的较重要因素；X5、X6、X7、X11和X12的解释力在0.01以下，虽一定程度影响土地利用的空间分异，但其单个因子对土地利用空间分异的影响较小。表明自然和经济因素均不同程度解释了保护区土地利用空间分异的驱动机制，其中X1、X2和X4等自然因素的解释力较高，故需加强保护区资源保护、生态保育和制度保障；同时，X8的解释力也较显著，应控制保护区居民点的规模，加强人类活动对生态系统破坏的风险管控。

交互探测器可探测驱动因子的交互作用对土地利用空间分异产生的影响。结果表明，任意两个因子的交互作用均大于单个因子的影响（图7），其中以X4与X8交互作用对土地利用空间分异的影响最为显著，解释力达0.56。交互作用对土地利用空间分异的解释力在0.40以上的还包括：X1∩X4（0.49）、X2∩X4（0.44）、X4∩X10（0.44）、X4∩X9（0.42）、X4∩X6（0.42）、X4∩X12（0.41）、X4∩X11（0.41）、X8∩X10（0.41）、X4∩X7（0.41）、X3∩X4（0.41）、X4∩X5（0.40）、X2∩X8（0.40）。其余因子交互类型的解释力均在0.400以下，但双因子交互较单因子对土地利用空间分异显示出更高的影响程度。从因子交互结果来看，X4与气象因子（X1和X2）、社会经济因子（X3、X5、X8和X9）和地形因子（X10、X11和X12）之间的交互作用明显强于自然因子内部交互和经济因子内部之间的交互作用。需依据保护区自然条件和社会经济发展水平，立足于保护自然的出发点，避免不合理人为土地利用干扰对区域生态系统的威胁。

3.4 土地利用的转型驱动因素

同样选取上述驱动因子，借助PLUS模型分析得到1990—2021年羌塘保护区土地利用变化的驱动因子贡献度。结果显示（图8），草地面积扩张的主导因子为X8，其贡献度为0.14，其次为X1（0.13）、X6（0.13）和X4（0.12）。从草地增加区域与X8的叠加图可看出（图9），草地增加区域集中于保护区北部和西部距居民点较远的地区，居民点较集中的南部地区草地增加量极少。水域面积扩张的主导因子为X6（0.26），其次为X4（0.17）、X10（0.12）和X2（0.01），增加区域集中于保护区西部、东北部和东南部湖泊连片带，距湖泊较远的高海拔地区水域面积增加较少。冰川积雪面积扩张的主导因子为X10（0.40），其次为X2（0.15）、X4（0.09）和X8（0.08），增加区域高度集中于西部昆仑山脉、中部唐古拉山脉和北部可可西里山脉等高寒地带。未利用地面积扩张的主导因子为X6（0.14），其次为X8（0.13）、X4（0.13）、X1（0.12）和X10（0.11），面积增加区域集中于保护区西部、东北部和东南部湖泊连片区，距湖泊较远的中部地区未利用地增加较少，可见季节交替导致的湖岸线变化及冬春季节的大风，使湖泊周边荒漠化现象加剧。1990—2021年羌塘保护区土地利用总体呈现草地、冰川积雪和未利用地面积缩减，水域面积不断扩张的趋势，其中草地、冰川积雪和未利用地转入水域的面积分别占水域增加面积的42.87%、4.50%和55.46%。因此，保护区水域面积的扩张在一定程度上解释了草地、冰川积雪和未利用地减少的原因。上述结果表明，相较于X8等经济因子，X1、X2、X4、X6和X10等自然因子对羌塘保护区土地利用扩张的影响更为显著。

进一步提取出主要土地利用类型面积扩张区域，并与高程数据相叠加（图10），发现驱动因子不仅影响着土地利用的空间变化，也制约着土地利用空间变化的形式。保护区西部、东北部和东南部湖泊集中的低洼地区，大面积草地和未利用地转为水域，促使湖泊面积不断扩张；中部唐古拉山脉南北麓的谷地，连片未利用地转为草地，未利用地收缩至山脚；北部地区湖泊周边有大量草地和水域转换为未利用地。未利用地转入冰川积雪的区域集中于西部昆仑山脉和唐古拉山脉等大型山脉高寒区。相较于南部居民点集中区，北部无人区内土地利用变化更为频繁，这进一步说明了羌塘保护区土地利用变化主要受制于自然因素，全球气候持续变暖的大背景下，需加强对保护区北部草地资源的保护，促进羌塘保护区生态系统的良性发展。

4 讨 论

4.1 土地利用的时空变化

根据土地利用动态变化分析可知，1990—2021年羌塘保护区土地利用主要呈现草地、未利用地、冰川积雪面积减少，水域面积不断扩张的趋势。其中，草地和未利用地转入水域的面积分别占水域增加面积的42.87%（2 384.36 km²）和55.46%（3 084.116 km²），同时水域增加面积集中于保护区东北部、东南部和西部的众多湖泊周边，即水域扩张主要通过湖泊淹没其周边未利用地和草地实现^[27]。从海拔和坡度来看，1990—2021年草地和未利用地的减少区域、水域增加区域均高度集中于海拔5 144 m、坡度16.09°的区域。在海拔5 144 m区间内，草地和未利用地减少面积分别占其减少总面积的86.94%（3 168.43 km²）和79.78%（1 978.12 km²），97.85%（5 448.88 km²）的水域增加面积集中于此。坡度16.09°的区域内，87.17%（2 491.37 km²）的草地和96.75%（1 948.55 km²）的未利用地减少面积也集中于此，91.05%（5 060.38 km²）的水域在此增加。此外，74.99%（1 057.73 km²）的冰川积雪减少面积集中于海拔5 262 m的区域，各级坡度内的减少量在88.80\sim191.30 km²之间。可见保护区土地利用变化主要发生在海拔5 144 m、海拔5 262 m和坡度16.09°的区域。从坡向来看，55.77%（1 594.09 km²）的草地减少面积集中于平地、东南、南和西南坡向，71.74%（936.67 km²）的未利用地减少于平面、东南、西和西北坡向，67.60%（953.70 km²）的冰川积雪减少于东北、东和西北坡向，水域除在平地增加较少外，其余各坡向均增加192.69～225.14 km²之间。可见阳坡上草地和阴坡上的未利用地和冰川积雪减少较显著，未来高海拔阴坡上的冰川积雪面临进一步退缩的趋势^[28]。

4.2 土地利用的驱动因子分析

根据羌塘保护区土地利用空间分异和面积变化的分析结果，可以发现气候、地形和距湖泊距离等自然因子为主要的影响因素。羌塘保护区为青藏高原腹地，人口稀疏，经济发展较为落后，居民点和公路极少。此外，高强度放牧虽影响草地生态系统稳定性^[7]，但羌塘地区放牧强度较小，对草场资源的干扰有限。故经济因子对羌塘保护区，乃至青藏高原主要土地利用类型变化的影响有限^[29]。同整个青藏高原一样^[29,30]，1990—2021年羌塘保护区周边4个主要气象站点（噶尔、申扎、班戈和那曲）的年均温明显上升，4地年均温由1990年的0.27℃上升至2021年的1.31℃。气候变暖影响下，青藏高原上的冰川一直在退缩^[31]，羌塘保护区西北部、东南部和东北部高大山脉的冰川积雪大量消融，融水注入湖泊，使众多湖泊得以扩张^[6]。可见全球气候变暖背景下，青藏高原地区脆弱的生态系统，极易受环境变化的影响^[32]。

从1990—2021年草地和未利用地的变化情况来看，草地流出面积主要转入水域和未利用地。保护区北部昆仑山脉南麓湖泊连片带周边，大面积草地退化为未利用地；保护区34°N以南的东南部湖泊聚集地带草地转为水域。前者是由于北部地区生态更为脆弱^[33]，随气温的不断上升，气候愈加干燥，大面积草场发生退化^[27]；后者是由于东南部地区植被覆盖度相对较高^[34]，在气温抬升下，高大山脉的冰川积雪加快融化，大面积草地资源被淹没。草地流入面积集中于中部唐古拉山脉、东北部可可西里山脉和西部昆仑山脉的山麓地带，主要由未利用地转入，得益于气温升高影响下的高大山脉冰雪融水的滋润，刺激了草地等植被的生长^[35]，局部生态环境有所改善。此外，保护区34°N以北的北部湖泊连片带，在四周高大山脉冰雪融水注入下，湖泊扩张淹没周围大片未利用地。

4.3 政策对土地利用的影响

为保护羌塘地区丰富的野生动物，1993年羌塘自然保护区由西藏自治区人民政府批准成立^[16]，促使1994—1995年草地面积增加0.56%（1 084.03 km²）、未利用地减少1.42%（960.60 km²）。2000年4月经国务院批准，羌塘保护区晋升为国家级保护区^[16]，使2001—2004年草地增长0.18%（340.03 km²）、未利用地缩减2.06%（1 444.70 km²）。2003年中国在内蒙古、西藏、青海等西部8省区全面实施退牧还草工程^[36]，后又于2011年完善退牧还草政策，促使2011—2021年羌塘保护区草地增长3.01%（5 648.52 km²）、未利用地减少9.78%（6 843.67 km²），与同时期在青藏高原东部和西北部的研究结果相似^[37]。设立自然保护区、实施退牧还草政策，通过禁牧封育、补播草种等方式，不仅使草原植被明显恢复^[29]，生态环境明显改善，也调动了广大农牧民保护和建设草原生态的积极性，增强了其草原保护意识。同时，退牧还草工程推行禁牧与休牧、舍饲与半舍饲相结合的生产方式，加快了传统草原畜牧业生产方式的转变。此外，在中央财政下拨补助、特色农牧产业与其他优势产业相结合、扶持牧区富余劳动力转产转业等举措推动下，退牧还草政策使农牧民收入稳步增长。设立保护区和退牧还草政策直接推动了羌塘保护区草地资源的恢复，对地区生态环境改善具有重要意义。

5 结 论

作为青藏高原的核心地带，羌塘保护区生态系统极易受环境变化影响。因此，监测土地利用变化对羌塘保护区生态稳定意义重大。利用1990—2021年土地利用数据，从地形角度分析了羌塘保护区土地利用的变化，并探测了自然和经济因子对土地利用变化和空间分异的影响，主要结论如下：

（1）羌塘保护区土地利用以草地为主导，未利用地、水域和冰川积雪也是较为重要的地类，其余地类分布极少。1990—2005年草地总体减少、未利用地呈增加趋势，2005—2021二者变化趋势相反。1990—2021年水域持续扩张，冰川积雪不断减少。草地、水域和未利用地三者转化频繁，水域增加部分主要由草地和未利用地转入。草地集中分布于保护区东南部、西部和东北部海拔4 939\sim5 405 m、坡度4.17°\sim25.93°的平坦地区；水域集中分布于保护区东南部、西部和东北部海拔4 939 m、坡度11.92°的平坦地区；冰川积雪集中于中部和西部海拔5 597 m、坡度16.09°的高大山脉，阴坡分布较多；未利用地主要分布于保护区北部和西部海拔4 818 m、坡度8.05°的区域，阳坡分布较多。

（2）草地在1990—2005年间仅于海拔5 889 m地区有所增加外，其余各级海拔、坡度和坡向均不同程度减少；2005—2021年间，除海拔4 939 m和海拔5 889 m的平地有所减少外，其余各尺度均有所增加。未利用地的变化趋势与草地正好相反，1990—2005年仅在海拔4 818 m和海拔5 889 m的平地区间有所减少，2005—2021年仅在坡度76.01°的区域有小幅增加。冰川积雪在2005—2021年除海拔5 889 m区域有少许增加外，1990—2021年其余各地形梯度均在减少。水域变化趋势与冰川积雪相反，1990—2021年仅在5 405 m海拔5 889 m的区域小幅减少，其余各地形梯度内均在增加。

（3）羌塘保护区地广人稀，经济活动对环境的影响有限，NDVI为土地利用空间分异的主导因子（q=0.40），年均温、年均降水和海拔等自然因子为重要因子。值得注意的是，与单因子相比，各驱动因子的交互对土地利用空间分异的解释力均显著增强。羌塘保护区主要土地利用类型的面积变化同样深受NDVI、距湖泊距离、海拔和年均温等自然因子影响，随着全球气候变暖，在冰川积雪的加速消融形成的冰雪融水不断注入下，保护区北部、东南部和西部的湖泊持续扩张，淹没了大量草地和未利用地。同时，北部地区干旱加剧，造成连片草地被未利用地侵占，荒漠化现象日益严重。

结 语

利用1990—2021年长时间序列逐年土地利用解译数据，不仅分析了羌塘保护区土地利用的动态变化，也从地形因子角度分析了土地利用在垂直梯度的动态变化，借助PLUS模型和地理探测器，探究了自然和经济因子对保护区土地利用面积变化及其空间分异的影响。土地利用地形梯度的动态变化，可为该地区制定因地制宜的生态环境保护策略提供参考。尽管如此，本研究仍有一定不足之处。首先，土地利用的动态变化分析高度依赖土地利用解译数据，不同分类标准和精度的解译数据会导致土地利用动态分析与实际情况产生偏差。采用总体精度为79.31%^[19]的土地利用数据，虽能进行可靠分析，但其误差不可避免，未来研究需进一步提高数据的准确性。其次，受数据获取源及数据精度限制，土地利用面积变化和空间分异的驱动因子仅选取12个，虽能进行可靠分析，但无法完全与实际情况相符，今后需进一步扩展影响因子，利用更高精度的气象数据来分析青藏高原地区土地利用变化尤为重要。最后，本文采用30 m分辨率的DEM数据提取坡度和坡向，与研究区实际的坡度和坡向存在一定失真。羌塘保护区土地利用动态变化分析结果表明，保护区需制定因地制宜的策略，如针对不同地形梯度内不同土地利用的变化特征，推出更为科学的保护措施。

参考文献

[1] 杜军,牛晓俊,袁雷,等.1971—2017年羌塘国家级自然保护区陆地生态环境变化[J].冰川冻土,2020,42(3):1017-1026.

[2] 杜军,高佳佳,次旺顿珠.1971—2019年羌塘国家级自然保护区日照时数变化特征及其影响因素[J].太阳能学报,2022,43(2):287-295.

[3] 吴彤,连新明,苗紫燕,等.藏羚行为的PAE编码系统[J].生态学杂志,2021,40(5):1449-1457.

[4] 魏子谦,徐增让,乔添,等.未来情景羌塘高原野牦牛栖息地分布及其影响因素[J]. 生态学报,2024,44(8):3317-3328.

[5] Hou Q Q, Ji Z X, Yang H, et al. Impacts of climate change and human activities on different degraded grassland based on NDVI[J]. Scientific Reports, 2022,12(1):15918.

[6] Liu J S, Wang S Y, Yu S M, et al. Climate warming and growth of high-elevation inland lakes on the Tibetan Plateau [J]. Global and Planetary Change, 2009,67(3):209-217.

[7] Lu X Y, Kelsey K C, Yan Y, et al. Effects of grazing on ecosystem structure and function of alpine grasslands in Qinghai-Tibetan Plateau: a synthesis[J]. Ecosphere,2017,8(1):e01656.

[8] 白淑英,陈灵梅,王莉,等.农牧交错区微地形对土地利用空间格局的影响研究[J].中国土地科学,2011,25(8):22-26+34+97.

[9] 何春阳,张金茜,刘志锋,等.1990—2018年土地利用/覆盖变化研究的特征和进展[J].地理学报,2021,76(11):2730-2748.

[10] 周建伟,罗君,马雪洋.拉萨河流域土地利用与生态系统服务价值时空演变及驱动因素[J].干旱区研究,2024,41(6):1021-1031.

[11] 周建伟,吴华,许童,等.西藏雅砻河流域生态系统服务价值地形梯度分异及驱动因素[J].水土保持研究,2025,32(2): 396-406.

[12] 冯彬,周建伟,于涛.冻融影响下生态系统服务价值和生态风险的时空演变特征及关联性——以南羌塘盆地东区为例[J]. 华南师范大学学报(自然科学版),2024,56(5):77-90.

[13] 胡书玲,余斌,卓蓉蓉,等.中国陆域地表人类活动与自然环境的空间关系研究[J].生态学报,2020,40(12):3935-3943.

[14] 陈运强,王荣远.土地利用空间格局与坡度坡向的关系[J].测绘与空间地理信息,2015,38(10):171-174+177.

[15] 杜军,周刊社,高佳佳,等.未来气候变暖对羌塘自然保护区高寒草地牧草青草期的影响[J].中国农业气象,2021,42(10): 845-858.

[16] Farrington J D, Tsering D. Snow leopard distribution in the Chang Tang region of Tibet, China[J]. GLOBAL Global Ecology and Conservation, 2020, 23: e01044.

[17] Feng F F, Yang Z S, Owens J R, et al. An assessment of endangered species habitat at large scale: chiru distribution across the Tibetan region of Chang Tang[J]. Folia Zoologica, 2016, 65(1): 65-71.

[18] Fox J L. Density of Tibetan antelope, Tibetan wild ass and Tibetan gazelle in relation to human presence across the Chang Tang Nature Reserve of Tibet, China[J]. Acta Zoologica Sinica, 2005, 51: 586-597.

[19] Yang J, Huang X. The 30 m annual land cover dataset and its dynamics in China from 1990 to 2019[J]. Earth System Science Data, 2021,13(8):3907-3925.

[20] 张学斌,罗君,石培基,等.格网尺度下张掖生态系统服务价值时空演变及地形梯度分异特征[J].应用生态学报,2020,31(2):543-553.

[21] 周建伟,吴华,赵鑫永,等.安宁河干旱河谷耕地地形梯度分异及其变化驱动因素[J].干旱地区农业研究,2024,42(4): 238-248.

[22] Qiu M L, Yuan C C, Yin G Y. Effect of terrain gradient on cadmium accumulation in soils[J]. Geoderma, 2020,375:114501.

[23] Xue L Q, Zhu B L, Wu Y P, et al. Dynamic projection of ecological risk in the Manas River basin based on terrain gradients[J]. Science of The Total Environment, 2019, 653(25): 283-293.

[24] Wang J F, Li X H, Christakos G, et al. Geographical Detectors-Based Health Risk Assessment and its Application in the Neural Tube Defects Study of the Heshun Region, China[J]. International Journal of Geographical Information Science, 2010, 24(1): 107-127.

[25] Liang X, Guan Q F, Clarke K C, et al. Understanding the drivers of sustainable land expansion using a patch-generating land use simulation (PLUS) model: A case study in Wuhan, China[J]. Computers, Environment and Urban Systems, 2021, 85: 101569.

[26] Yao Y, Liu X P, Li X, et al. Simulating urban land-use changes at a large scale by integrating dynamic land parcel subdivision and vector-based cellular automata[J]. International Journal of Geographical Information Science, 2017, 31(12): 2452-2479.

[27] Cui X F, Graf H. Recent land cover changes on the Tibetan Plateau: a review[J]. Climatic Change, 2009, 94(1): 47-61.

[28] Xia M, Jia K, Zhao W W, et al. Spatio-temporal changes of ecological vulnerability across the Qinghai-Tibetan Plateau[J]. Ecological Indicators, 2021, 123: 107274.

[29] Wang S S, Tan X, Fan F L. Landscape Ecological Risk Assessment and Impact Factor Analysis of the Qinghai-Tibetan Plateau[J]. Remote Sensing, 2022, 14(19): 4726.

[30] Chen H, Zhu Q A, Peng C H, et al. The impacts of climate change and human activities on biogeochemical cycles on the Qinghai-Tibetan Plateau[J]. Global Change Biology, 2013, 19(10): 2940-2955.

[31] You Q L, Wu T, Shen L C, et al. Review of snow cover variation over the Tibetan Plateau and its influence on the broad climate system[J]. Earth-Science Reviews, 2020, 201: 103043.

[32] Qiu P F, Wu N, Luo P, et al. Analysis of dynamics and driving factors of wetland landscape in Zoige, Eastern Qinghai-Tibetan Plateau[J]. Journal of Mountain Science, 2009, 6(1): 42-55.

[33] Liang S, Chen J, Jin X M, et al. Regularity of vegetation coverage changes in the Tibetan Plateau over the last 21 years[J]. Advances in Earth Science, 2007, 22(1): 33-40.

[34] Zhang L N, Zhang H Q, Xu E Q. Information entropy and elasticity analysis of the land use structure change influencing eco-environmental quality in Qinghai-Tibet Plateau from 1990 to 2015[J]. Environmental Science And Pollution Research, 2022, 29(13): 18348-18364.

[35] Huang K, Zhang Y J, Zhu J T, et al. The influences of climate change and human activities on vegetation dynamics in the Qinghai-Tibet Plateau[J]. Remote Sensing, 2016, 8(10): 876.

[36] Liu X, Ma Z W, Huang X T, et al. How does grazing exclusion influence plant productivity and community structure in alpine grasslands of the Qinghai-Tibetan Plateau? [J]. Global Ecology Conservation, 2020, 23:e01066.

[37] Diao C, Liu Y, Zhao L, et al. Regional-scale vegetation-climate interactions on the Qinghai-Tibet Plateau[J]. Ecological Indicators, 2021, 65:101413.