中国典型山地土地利用变化与资源生态效应*

戴尔阜，王亚慧，马良，李双成，张红旗，辛良杰，徐尔琪，高江波， 朱连奇，王玉宽

①中国科学院地理科学与资源研究所，北京 100101；②中国科学院大学，北京 100049；③中国科学院陆地表层格局与模拟重点实验室，北京 100101；④北京大学 城市与环境学院 地表过程分析与模拟教育部重点实验室，北京 100871；⑤河南大学 环境与规划学院，河南 开封 475004；⑥中国科学院水利部成都山地灾害与环境研究所，成都 610041

山区以其独特的气候、水文、土壤、地质条件孕育着森林、草地、湿地等不同类型的生态系统，其往往是大江大河的发源地及上游地区，所发挥的调节气候、水源涵养、净化空气等生态功能不仅惠及山区居民，同时也是山区河流下游区域以及更为广阔区域居民福祉的重要来源[1]。同时，山区因其陡峭的地形使其物流、能流具有不稳定性，在气候变化以及人类活动强度增加的背景下，山区生态系统将更为脆弱[2]。在中国，山区面积广大，约占到陆域国土面积的70 %，同时有45 %的人口生活在山区，人口密度较高[3]。中国山地，特别是西部山地，地形复杂、海拔梯度变化大、耕地分散、土壤层浅薄、土地利用空间差异较大，同时在山区城镇化建设、生态工程、生态移民、灾害避险搬迁等驱动力的影响下，山区土地利用的方式及强度发生了剧烈的变化，并通过水文、地球化学、生物等过程变化进而影响生态系统服务的发挥以及人类福祉的获得。在这种山区空间异质性背景下，关注资源生态效应及其变化对我国山区可持续发展具有重大战略意义。本研究分析了中国太行山区、横断山区以及黔桂喀斯特山区三个典型山区土地要素的空间变异特征及其功能，甄别了不同山地关键生态系统服务类型，并阐明了其资源与生态效应。

1 中国典型山地土地利用结构变化

在土地利用结构上，三大山区表现为一致性的特征为耕地、林地、草地三者所占比例均达到90 %以上(图1)，但在不同土地利用类型面积比例的排序上有所不同。其中太行山区位于我国第二阶梯东沿平原到山地的过渡带，是典型的土石山区，土壤贫瘠[4]，历史上人类活动对天然植被造成了严重破坏[5]。太行山区三种土地利用类型面积排序为耕地＞林地＞草地，人均耕地为0.10 hm2/人，人均林地为0.09 hm2/人，人均草地为0.08 hm2/人。横断山区位于青藏高原东南缘，横跨我国一二三级地形阶梯，是我国西南林区的主要组成部分，林地占比较大[6-7]。横断山区三种土地利用类型面积排序为林地＞草地＞耕地，人均林地为1.47 hm2/人，人均草地为1.05 hm2/人，人均耕地为0.26 hm2/人，耕地资源短缺。黔桂喀斯特山区多以石山形式出现，又称石山地区、石灰岩区或岩溶区，土地破碎程度较高[8]，其三种土地利用类型面积排序为林地＞耕地＞草地，人均林地面积为1.22 hm2/人，人均耕地为0.56 hm2/人，人均草地为0.31 hm2/人。从时间变化来看，三大山区建设用地表现为持续增加，耕地为先增后减，而林地、草地则表现为先减后增的趋势。

图1 典型山区土地利用(左图：太行山区；中图：横断山区；右图：黔桂喀斯特山区)

表1 典型山区人均土地资源 hm2/人

2 典型山地土地利用垂直及水平空间异质性

三大山区的土地利用情况具有明显的垂直和水平空间异质性。在垂直方向上，受海拔高度影响，导致气候、土壤以及植被分布出现垂直差异，从而形成土地利用/覆被垂直空间异质性。在水平方向上，人口的增加、城市化进程的加快以及频繁的社会经济活动对山区城市扩张有显著影响，随着到城市中心距离的增加，人类活动影响下的土地利用具有一定的水平空间异质性。

2.1 土地利用垂直空间异质性

垂直方向上，将高程按照等间距进行区域划分并统计主要土地利用类型(耕地、林地、草地)所占比例，结果表明该比例随着海拔高度的增加呈现明显的规律性变化(图2)。横断山区：耕地集中分布于3 500 m以下，随着海拔的增加，比例逐渐减小；林地主要集中分布于4 500 m以下，林地面积随着海拔高度增加呈现先增加后减少的趋势，在3 400 m时达到峰值，为单峰形态；草地在各海拔高度均有分布。在时间变化上，横断山区在2000—2010年间耕地面积明显较少，主要位于1 300 m高程以下地区；林地在此期间增加的区域与草地面积减少的区域在垂直方向上具有一致性，主要位于4 000 m左右。黔桂喀斯特山区：耕地在低海拔分布的比例最高，在海拔800～1 400 m分布比例出现第二个峰值；林地在低于1 200 m海拔的各等级分布较大，大于1 200 m海拔之后迅速减少；草地在垂直方向上分布则更加均匀，1 200～1 400 m等级的草地面积比例最大。在时间变化上，黔桂喀斯特山区 0～200 m海拔耕地减少面积远大于增加面积，600～1 400 m耕地增加和减少面积出现峰值；草地和林地在1 000～1 400 m出现变化的峰值，两类用地转化分布相对应。

图2 典型山区不同海拔高度耕地、林地、草地面积比例

2.2 土地利用水平空间异质性

土地空间异质性不仅包括土地的组成，同时也包含着土地结构的异质性[9-10]。以横断山区为例，通过对不同城市建立“城-乡”梯度，并结合景观指数分析各缓冲区内土地结构的变化情况，结果表明：从城市中心到外围，在各个缓冲区内，城市用地逐渐被其他用地类型代替，缓冲区内景观结构发生相应的变化，景观LPI指数和CONTAG指数呈现减小的趋势，而且在城市用地与其他用地交错程度最高的缓冲区内，景观优势度与连通度较低，LPI指数和CONTAG指数达到低值，之后两者呈现增加趋势；而AWMFD指数和SHDI指数随着到城市中心距离的增加先呈现增加的趋势，而在城市用地与其他用地交错程度最高的缓冲区内达到高值，之后两者呈减小趋势(图3)。

图3 缓冲区内景观指数随到城市中心距离的变化

3 中国典型山区资源生态效应

生态系统服务(ecosystem service)是人类从生态系统中获得的各种利益，包括产品供给与服务提供两方面[11-12]。千年生态系统评估(MA)将其概括为供给服务、调节服务、文化服务和支持服务4种类型[13]，同时也是资源生态效应的重要指标。结合山区垂直带主导资源提供及生态系统服务类型，采用InVEST(综合服务价值权衡工具)模型对典型山区产水量、土壤保持、碳储存和生境质量4种关键生态系统服务进行量化评估。

(1)产水量

基于 Budyko 水热耦合平衡假设计算生态系统总的产水量(water yield)[14-15]。根据流域水量平衡原理，流域内多年平均降水量是径流量、蒸散发量和流域蓄水变化量的和。对于一个自然闭合流域，在没有显著的人类活动影响下，多年平均需水量趋于平衡。因此，流域产水量可以看作降水量与蒸发量的差值[16]。三大典型山区产水量评估结果如图4所示。2010年，太行山区中部蒸散发量大于降水量，产水量为负值，南北两侧降水量大于蒸散发量，年均产水量是75 mm/a；横断山区的年均产水量约为385 mm/a，产水量具有明显的空间差异，与降雨量的分布格局大致相同，研究区西南部和东部西川盆地边缘区域为高值区；喀斯特山区产水量在262.78～1 707.45 mm之间，年均产水量为968.10 mm/a，中部和东南地区产水量较大，北部地区产水量相对较小。

(2)土壤侵蚀与保持

图4 三大典型山区2010年产水量变化(左图：太行山区；中图：横断山区；右图：黔桂喀斯特山区)

通用土壤流失方程(RUSLE)综合考虑影响土壤侵蚀的降雨侵蚀力、土壤可蚀性、地形、植被覆盖以及管理因子的影响而得到了广泛的应用[17]。InVEST 模型在其基础上考虑了泥沙输送率(sediment delivery ratio, SDR)使得计算结果更符合实际[18]，而潜在土壤侵蚀(不考虑植被覆盖)与实际土壤侵蚀量之间的差值即可以看作土壤保持量。三大典型山区土壤侵蚀/保持量评估结果如图5所示。2010年，太行山区土壤侵蚀年平均值为65.5 t·km-2·a-1，土壤侵蚀量主要受到地形因子的影响，坡度较大的地区侵蚀量较大；横断山区的土壤侵蚀年平均值约为2 154.3 t·km-2·a-1，土壤保持量年平均值为3 897.6 t·km-2·a-1，横断山区的南部地区土壤保持量略高于北部；黔桂喀斯特山区土壤保持量年平均值为1 035.90 t/hm2，相对来说，中部和西北地区土壤保持量较大，而北部、南部和东南地区土壤保持量较小。

(3)碳储量

碳储量综合考虑了地上生物量、地下生物量、土壤和枯枝落叶四大碳库。横断山区和黔桂喀斯特山区碳储量评估结果如图6所示。2010年，横断山区的碳储量均值为29.5 t/hm2。碳储量的空间分布受土地利用格局的影响存在着显著的空间差异，横断山区南部地区碳储量整体高于北部地区；从土地利用类型上来看，林地的碳储量最高，其他土地利用类型的碳储量相对较低。黔桂喀斯特山区的平均碳储量为30.75 t/hm2，西北部和东南部呈现较低的分布，中部较高。

图5 三大典型山区2010年土壤保持(左图：太行山区；中图：横断山区；右图：黔桂喀斯特山区)

图6 横断山区(左图)和黔桂喀斯特山区(右图)2010年碳储量

(4)生境质量

生境是指一个特定生物体占有(包括生存和繁殖)的地区的资源和条件[19]。生境质量指的是生态系统提供适合个体和种群维持的条件，在模型中作为连续变量考虑，基于可利用的生存、繁殖和种群维持所需的资源，从低到高划分[19]。通常来说，生境质量降低反映了附近区域土地利用强度的增加[20-22]。生境质量是评估生物多样性的一个重要指标。这里的生境质量主要指的是胁迫源对生境的影响程度。生境质量由4个因素决定：①不同胁迫源的相对影响；②不同生境对胁迫源的敏感性；③生境与胁迫源的距离；④生境被保护的程度[16]。

横断山区和黔桂喀斯特山区生境质量评估结果如图7所示。生境质量是山区一项重要的生态服务功能，区域生境质量反映了区域的生物多样性水平。2010年，横断山区生境质量指数为0.653 1。中西部地区生境质量较高；东南部生境质量较差，主要是受耕地、城镇及人类活动影响；北部地区，主要是受公路及道路的影响，为生境质量较弱的现状分布区。黔桂喀斯特山区的平均生境质量指数为0.571 3，高值区位于中部地区。

4 30年来中国山地地表绿色覆被变化及其对气候因子的响应

植被是连接土壤、大气和水分等的纽带，在调节水量、保持土壤、减缓温室气体浓度上升、维持生物多样性及整个生态系统稳定等方面都具有十分重要的作用[23]。众多反映植被状况的因子中，归一化植被指数(NDVI)是监测地区或全球植被与生态环境的有效指标，能很好地反映植被覆盖和生物量等的变化[24-25]。NDVI 与植被覆盖度、叶面积指数、生物量和生产力等性状关系密切。20世纪90年代以来，全球环境发生了剧烈的变化，NDVI动态变化与气候因子的关系成为了地理学和生态学的研究热点[26-28]。

图7 横断山区(左图)和黔桂喀斯特山区(右图)2010年生境质量

采用Theil-Sen斜率计算方法[29-31]，通过提取中国山地并对1982—2010年山地绿色覆被变化进行分析，得到中国山地1982—2012年均NDVI变化率为0.000 32，整体呈现“绿化”态势(图8)。其中秦巴山地以北至燕山山地、云贵高原和台湾山地，“绿化”趋势显著，东北山地“棕化”态势明显。山地地表覆被变化程度上，绝大部分区域为稳定少变，其次是显著改善，其他依次为轻微改善、显著退化和轻微退化。以1998年作为拐点，中国山地大部分地区植被变绿趋势呈现减缓的趋势[32-33]。

图8 1982—2012年中国山地地表覆被变化趋势及变化比例

综合考虑山地地表绿色覆被变化对气候变化影响的滞后效应基础上，计算植被生长与温度、降水和太阳辐射的偏相关关系，结果表明中国山地生长季NDVI与温度因子、降水因子以及辐射因子的平均偏相关系数分别为0.347、0.136和0.250，说明山区植被受到温度和辐射的影响大于降水的影响(图9)。

5 结论

山区是人地关系最为复杂的区域之一，人口密度大而可利用土地有限是山区面临的重要问题，协调山区经济发展及生态屏障功能是山区当前发展所面临的巨大挑战。土地利用作为人类活动的集中体现，在促进经济水平发展的同时也影响着生态功能的发挥，量化土地空间异质性以及相应生态功能是维持山区平衡发展的关键。本文以三大山区为对象分析了其土地组成及结构的空间变化，结果表明：三大山区在土地利用组成方面均以耕地、林地、草地为主导类型；在土地利用结构方面表现为垂直与水平空间异质性并存的特征。垂直方向上，主要土地利用类型占比随海拔的增高呈现出规律性变化；而在水平方向上，景观结构表现为随到城市中心距离的增加表现为梯度变化，同时城乡交错地带是景观最为复杂的区域。在此基础上，甄别了三大山区的主导生态系统服务(产水量、土壤保持、碳储量、生境质量)并进行了量化评估。最后，针对中国山地地表覆被变化的研究表明：近30年来中国山地整体呈现“绿化”态势，其中1998年之后植被变绿趋势减慢；在考虑时滞效应基础上，中国山地植被受到温度和辐射的影响大于降水的影响。本文的研究结果对合理制定土地规划利用政策具有重要意义。未来通过对山区不同情景下的土地及生态系统服务进行分析，量化各利益相关方的经济收支平衡，将能为制定生态补偿标准以及山区经济与生态协调发展提供更为具体的科学依据。

图9 中国山地植被生长季NDVI与气候因子的偏相关系数空间分布图(左图：NDVI与温度；中图：NDVI与太阳辐射；右图：NDVI与降水)

(2018年1月3日收稿)

[1] LOCATELLI B, LAVOREL S, SLOAN S, et al. Characteristic trajectories of ecosystem services in mountains[J]. Frontiers in Ecology and the Environment, 2017, 15(3). doi: 10.1002/fee.1470.

[2] 钟祥浩. 中国山地生态安全屏障保护与建设[J]. 山地学报, 2008,26(1): 2-11.

[3] 邓伟, 唐伟. 试论中国山区城镇化方向选择及对策[J]. 山地学报,2013, 31(2): 168-173.

[4] 刘欣, 葛京凤, 冯现辉. 河北太行山区土地资源生态安全研究[J]. 干旱区资源与环境, 2007, 21(5): 68-74.

[5] 李晓荣, 高会, 韩立朴, 等. 太行山区植被NPP时空变化特征及其驱动力分析[J]. 中国生态农业学报, 2017, 25(4): 498-508.

[6] 沈文娟, 徐婷, 李明诗. 中国三大林区森林破碎化及干扰模式变动分析[J]. 南京林业大学学报(自然科学版), 2013, 37(4): 75-79.

[7] 王强, 张廷斌, 易桂花, 等. 横断山区2004—2014年植被NPP时空变化及其驱动因子[J]. 生态学报, 2017, 37(9): 3084-3095.

[8] 徐卓, 刘芸, 尹小华, 等. 西南喀斯特山地林下经济发展模式研究[J].西南师范大学学报(自然科学版), 2012, 37(11): 59-65.

[9] FAHRIG L, NUTTLE W K. Population ecology in spatially heterogeneous environments [M]//Ecosystem Function in Heterogeneous Landscapes. New York: Springer, 2005: 95-118.

[10] JONES K B, ZURLINI G, KIENAST F, et al. Informing landscape planning and design for sustaining ecosystem services from existing spatial patterns and knowledge [J]. Landscape Ecology, 2013, 28(6):1175-1192.

[11] COSTANZA R, D'ARGE R, DE GROOT R, et al. The value of the world's ecosystem services and natural capital [J]. Nature, 1997,387(6630): 253-260.

[12] DAILY G. Nature's services: societal dependence on natural ecosystems [M]. Washington, D C: Island Press, 1997.

[13] ASSESSMENT M E. Ecosystem and human well-being: general synthesis [M]. Washington, D C: World Resource Institute, 2005.

[14] BUDYKO M. Evaporation under natural conditions [M]. English translation by IPST, Jerusalem. Leningrad: Gidrometeorizdat, 1948.

[15] 傅抱璞. 论陆面蒸发的计算[J]. 大气科学, 1981, 5(1): 23-31.

[16] SHARP R, TALLIS H T, RICKETTS T, et al. InVEST 3.2.0 User’s Guide [M]. The Natural Capital Project, Stanford University,University of Minnesota, The Nature Conservancy, and World Wildlife Fund, 2015.

[17] RENARD K G, FOSTER G R, WEESIES G A, et al. Predicting soil erosion by water: a guide to conservation planning with the revised universal soil loss equation (RUSLE) [M]// USDA Agriculture Handbook No. 703, 1997.

[18] VIGIAK O, BORSELLI L, NEWHAM L, et al. Comparison of conceptual landscape metrics to define hillslope-scale sediment delivery ratio [J]. Geomorphology, 2012, 138(1): 74-88.

[19] HALL L S, KRAUSMAN P R, MORRISON M L. The habitat concept and a plea for standard terminology [J]. Wildlife Society Bulletin,1997, 25(1): 173-182.

[20] NELLEMANN C, KULLERUD L, VISTNES I, et al. GLOBIO.Global methodology for mapping human impacts on the biosphere,UNEP/DEWA/TR. 01-3 [R]. United Nations Environment Programme,2001.

[21] MCKINNEY M L. Urbanization, biodiversity, and conservation [J].BioScience, 2002, 52(10): 883-890.

[22] FORMAN R T. Road ecology: science and solutions [M]. Washington,D C: Island Press, 2003.

[23] 孙红雨, 王长耀, 牛铮, 等. 中国地表植被覆盖变化及其与气候因子关系──基于NOAA时间序列数据分析[J]. 遥感学报, 1998, 2(3):204-210.

[24] 范娜, 谢高地, 张昌顺, 等. 2001年至2010年澜沧江流域植被覆盖动态变化分析[J]. 资源科学, 2012, 34(7): 1222-1231.

[25] 吕洋, 董国涛, 杨胜天, 等. 雅鲁藏布江流域NDVI时空变化及其与降水和高程的关系[J]. 资源科学, 2014, 36(3): 603-611.

[26] BRUNSDON C, FOTHERINGHAM S, CHARLTON M.Geographically weighted regression [J]. Journal of the Royal Statistical Society: Series D (The Statistician), 1998, 47(3): 431-443.

[27] WALTHER G R, POST E, CONVEY P, et al. Ecological responses to recent climate change [J]. Nature, 2002, 416(6879): 389-395.

[28] LESAGE J P. A family of geographically weighted regression models[M]//Advances in Spatial Econometrics. Berlin: Springer, 2004: 241-264.

[29] SEN P K. Estimates of the regression coefficient based on Kendall's tau [J]. Journal of the American Statistical Association, 1968, 63(324):1379-1389.

[30] THEIL H. A rank-invariant method of linear and polynomial regression analysis [M]//RAJ B, KOERTS J (eds). Henri Theil’s Contributions to Economics and Econometrics. Advanced Studies in Theoretical and Applied Econometrics, Vol. 23. Dordrecht: Springer, 1992: 345-381.

[31] 匡薇, 马勇刚, 李宏, 等. 中亚1999—2012年间土地退化强度与趋势分析[J]. 国土资源遥感, 2014, 26(4): 163-169.

[32] LIU Y, LI S. Spatial and temporal patterns of global NDVI trends:correlations with climate and human factors [J]. Remote Sensing,2015, 7(10): 13233-13250.

[33] 韩雅, 朱文博, 李双成. 基于GWR模型的中国NDVI与气候因子的相关分析[J]. 北京大学学报(自然科学版), 2016, 52(6):1125-1133.