人类的足迹遍布整个星球,地球进入了人类世(Kareiva et al., 2007; Steffen et al., 2007)。然而,地球的承载能力是有限的。随着对气候变化和环境危机理解的加深,人类与自然之间的矛盾被认为是阻碍可持续发展的根本原因(Anderson et al., 2019; Dorninger et al., 2017)。人类活动在生物圈完整性以及氮和磷的循环方面已超出了行星边界,生活在一个安全和公正空间的挑战依然巨大(O’Neill et al., 2018)。因此,理解社会经济发展与生态环境演变之间的关系对于全球可持续发展至关重要。此外,平衡人类活动与生态环境是许多联合国倡议的核心。
写作 - 原始草稿:Haimeng Liu, Yi Cheng
写作 - 审阅与编辑:刘海萌,刘志峰,李启瑞,张海燕
可持续发展目标(SDGs),包括目标 6、7、11、12、13、14 和 15(Cheng et al., 2021;Liu, Hull, et al., 2018),以及许多全球组织或项目也探讨了人类-自然系统的耦合机制和协同发展方法(Kates et al., 2001;Liu, Fang, & Fang, 2020)。分析人类活动与生态环境之间的矛盾与协调对人类的命运和福祉具有重要的现实意义。
青藏高原(QTP)平均海拔超过 4000 米,广泛分布着冰川和永久冻土,是亚洲九大河流的源头,为超过 15 亿人提供淡水、食物和其他生态系统服务,被称为地球的第三极和亚洲水塔(Yao et al., 2012)(图 1)。它在全球气候调节和生物多样性维持中发挥着重要作用(Li et al., 2020)。然而,由于其高海拔、寒冷的气温和干燥的环境,QTP 的生态系统对人类活动极为敏感和脆弱(Han et al., 2022;Liu, Milne, et al., 2018)。从 1980 年到 2020 年,人口增加了 1.56 倍,城市建设面积增加了近三倍,高速公路长度增加了四倍以上(Zhang et al., 2019)。过度放牧等强烈的人类活动对当地草原生态系统施加了巨大的压力,这种压力可能远远超过气候变化的影响(Wei et al., 2022)。人类活动的加剧对 QTP 生态系统的稳定性和可持续发展目标(SDGs)的实现构成了重大威胁(Yang et al., 2022)。 因此,阐明 QTP 上的人类与自然的关系对地方和全球可持续发展具有重要意义。
人类与自然之间的互动自人类首次出现以来就存在(Costanza et al., 2007)。在环境问题日益严重的背景下,许多来自不同学科背景的学者专注于这一与可持续发展密切相关的领域,并产生了宝贵的见解(Harden, 2012)。研究主题包括人类活动对自然的正面和负面影响(Kareiva et al., 2007)、自然资源和生态系统对社会发展的约束和承载能力(Fan et al., 2017),以及它们之间的互动、反馈、权衡和适应(Comberti et al., 2015;Holling & Gunderson, 2002;Kareiva et al., 2007)。一些研究表明,减少人类活动的规模对于缓解社会危机是必要的(Dirzo et al., 2022);另一些研究则表明,如果我们在生产方式、能源使用、生态保护等方面进行适当的改变,可以实现双赢的愿景(Cao et al., 2009;Tallis et al., 2018)。
在这个研究领域,城市化地区(Fang et al., 2016)、流域(Yin et al., 2021)、森林(Davidson et al., 2012)和沿海地区(Yi et al., 2018)受到了更多关注。对于自然系统,研究人员关注生物多样性(Cardinale et al., 2012)、植被(Zhang, Yang, et al., 2022)、空气污染(Liu, Cui, & Zhang, 2022)、河流环境(Leprieur et al., 2008)、土壤侵蚀(Wei et al., 2006)等。土地利用和覆盖变化(Cheng et al., 2022)、人口密度(Cropper & Griffiths, 1994)和夜间光照(Cai et al., 2021)通常被用作人类活动的代理变量;一些学者还构建了一个名为人类足迹的综合指数(Li, Wu, Gong, et al., 2018;Venter et al., 2016)。用于分析人类与自然之间复杂关系的方法包括网络方法(Kluger et al., 2020)、系统动力学(Reed et al., 2022;Yoon et al., 2021)、耦合协调分析(Jiang et al., 2022;Shi, Feng, et al., 2022)、复杂系统理论(Bennett & McGinnis, 2008;Wang & Grant, 2021)、环境库兹涅茨曲线模型(Zhou et al.,2015),元耦合框架(Liu,2017),供需平衡分析(Liu,Xing 等,2022),问卷调查(Allison 等,2021;Liu,Zhang 等,2020),以及定性分析(Ives 等,2017)。这些方法从不同角度推动了复杂人类-自然系统的研究,但在细网格尺度上面临重大限制。某些使用栅格数据的研究:Vačkár 等(2012)使用 10 公里网格分析了捷克共和国人口密度、土地利用强度和生物多样性之间的空间关系。Shi,Feng 等(2022)利用耦合协调度阐明了中国 10kmxx10km10 \mathrm{~km} \times 10 \mathrm{~km} 地区景观生态风险与城市化之间的协调关系。Zhou 等(2023)测量了中国城市中城市化与植被覆盖之间的冲突与协调动态,使用了 250 米的植被指数。Simkin 等(2022)预测了城市土地扩张对三种共享社会经济路径(SSP)情景下生物多样性的未来影响。然而,这些研究集中于某些人类活动或自然元素,缺乏对人类和自然系统的整体考虑。
古老的人类自 16 万年前起就占据了青藏高原(Chen et al., 2019)。在过去的半个世纪里,人类活动对青藏高原生态系统造成了有限但迅速增加的干扰,特别是在草地、宝贵的水源保护区和生物多样性保护区的最大比例上(Jie et al., 2015;Li, Zhang, Wang, et al., 2018)。关于青藏高原上人类与自然关系的研究
图 1. 青藏高原的自然地理。
主要关注城市化、土地利用、基础设施建设、农业和畜牧业对生态环境的影响(Miao et al., 2021; Wang et al., 2015; Zhang, Zhang, et al., 2022)。这里的城市化水平相对较低,但近年来城市人口迅速增长,建设用地快速扩张(Wang et al., 2020)。尽管城市化会导致更多的污染物排放,但通过将人口从农业和牧区转移,可以减少人类活动对生态系统造成的整体干扰(Tian et al., 2021)。2000 年至 2015 年间,城市化与生态环境之间的协调性呈上升趋势(Feng & Li, 2021)。草地退化主要是由于短期内人类活动造成的,这进一步导致生物多样性丧失、土壤侵蚀加剧、碳汇减少以及当地人类福祉下降(Chen et al., 2020; Dong et al., 2020)。青藏铁路和川藏铁路的建设和运营将给当地生态环境带来潜在风险(Cui et al., 2022; Luo et al., 2020; Qin & Zheng, 2010)。 从系统的角度,利用承载能力模型(Niu et al., 2020)分析了产业、人口与资源环境的共同发展,通过将人类活动与生态系统服务联系起来评估生态质量(Sun et al., 2020;Wang, Huang, et al., 2022),并通过整合生态能力、生态足迹和资源效率来衡量区域可持续性(Fan & Fang, 2022)。与此同时,中国在青藏高原的生态保护与恢复项目,包括植树造林、草原恢复和设立自然保护区,正在改善许多地区的生态环境(Fu et al., 2021;Jin et al., 2019)。近年来,矿业和城市的污水排放和废气排放、因游客增加而产生的垃圾以及光污染已成为新的严重环境问题(Chen, Zhang, et al., 2021;Feng et al., 2020;Wang, Lv, et al., 2022)。
存在许多度量标准来量化人类活动的强度(人类足迹,人类改造),这些标准常用于描述人类对生态系统的压力(Venter et al., 2016; Williams et al., 2020)。本研究构建的 HAI 在很大程度上遵循了先前研究建立的评估指标体系,包括人口密度、GDP、夜间灯光、土地覆盖和道路网络(Huang et al., 2022; Mu et al., 2022; Xu et al., 2016)。然而,考虑到草原占据了 60%60 \% 的 QTP(Li
Assigned from 0 to 100 according to intervals determined by 10 equal quantiles
(Mu et al., 2022)| Assigned from 0 to 100 according to intervals determined by 10 equal quantiles |
| :--- |
| (Mu et al., 2022) |
Data set Year and format Data source Data processing
Population density " 2000, 2005, 2010, 2015, 2020;
Raster, 1km" WorldPop For all locations with more than 1,000 people km^(2), we assigned a score of 100 . For more sparsely populated areas, we logarithmically scaled the score using 33.333 xx log (population density +1 ) (Williams et al., 2020)
GDP density " 2000, 2005, 2010, 2015, 2019;
Raster, 1km" Chen et al. (2022) Scale data into the 0-100 range using Min-Max Normalization, and the maximum value is the average of the top 50 values in the 5 years
Nighttime light " 2000, 2005, 2010, 2015, 2020;
Raster, 1km" Zhang, Ren, et al. (2021) "Assigned from 0 to 100 according to intervals determined by 10 equal quantiles
(Mu et al., 2022)"
Grazing intensity " 2000, 2005, 2010, 2015, 2019;
Raster, 250m" Liu (2021) The data in 2005 used the average of 2000 and 2010, and the data in 2015 used the average of 2010 and 2019. Resampling to 1-km resolution. Scale data into the 0-100 range using Min-Max Normalization, and the maximum value is the average of the top 50 values in the 5 years
Land use/cover " 2000, 2005, 2010, 2015, 2020;
Raster, 1km" www.resdc.cn/data. aspx?DATAID=335 Assign score: 100 for urban land, 50 in 1 km buffer; 80 for construction land, 40 in 1 km buffer; 60 for cropland, 20 in 1 km buffer; 0 for other lands
Road network " 2000, 2005, 2010, 2015, 2018;
Shapefile " Digitization by authors Assign a score of 60 for 0.5 km out for either side of a highway, national road, and provincial road. Assign 80 and 40 for 0.5 and 1.5 km out for either side of a railway (Venter et al., 2016). Convert to raster data| Data set | Year and format | Data source | Data processing |
| :---: | :---: | :---: | :---: |
| Population density | $\begin{aligned} & \text { 2000, 2005, 2010, 2015, 2020; } \\ & \text { Raster, } 1 \mathrm{~km} \end{aligned}$ | WorldPop | For all locations with more than 1,000 people $\mathrm{km}^{2}$, we assigned a score of 100 . For more sparsely populated areas, we logarithmically scaled the score using $33.333 \times \log$ (population density +1 ) (Williams et al., 2020) |
| GDP density | $\begin{aligned} & \text { 2000, 2005, 2010, 2015, 2019; } \\ & \text { Raster, } 1 \mathrm{~km} \end{aligned}$ | Chen et al. (2022) | Scale data into the $0-100$ range using Min-Max Normalization, and the maximum value is the average of the top 50 values in the 5 years |
| Nighttime light | $\begin{aligned} & \text { 2000, 2005, 2010, 2015, 2020; } \\ & \text { Raster, } 1 \mathrm{~km} \end{aligned}$ | Zhang, Ren, et al. (2021) | Assigned from 0 to 100 according to intervals determined by 10 equal quantiles <br> (Mu et al., 2022) |
| Grazing intensity | $\begin{aligned} & \text { 2000, 2005, 2010, 2015, 2019; } \\ & \text { Raster, } 250 \mathrm{~m} \end{aligned}$ | Liu (2021) | The data in 2005 used the average of 2000 and 2010, and the data in 2015 used the average of 2010 and 2019. Resampling to $1-\mathrm{km}$ resolution. Scale data into the $0-100$ range using Min-Max Normalization, and the maximum value is the average of the top 50 values in the 5 years |
| Land use/cover | $\begin{aligned} & \text { 2000, 2005, 2010, 2015, 2020; } \\ & \text { Raster, } 1 \mathrm{~km} \end{aligned}$ | www.resdc.cn/data. aspx?DATAID=335 | Assign score: 100 for urban land, 50 in 1 km buffer; 80 for construction land, 40 in 1 km buffer; 60 for cropland, 20 in 1 km buffer; 0 for other lands |
| Road network | $\begin{aligned} & \text { 2000, 2005, 2010, 2015, 2018; } \\ & \text { Shapefile } \end{aligned}$ | Digitization by authors | Assign a score of 60 for 0.5 km out for either side of a highway, national road, and provincial road. Assign 80 and 40 for 0.5 and 1.5 km out for either side of a railway (Venter et al., 2016). Convert to raster data |
rho1\rho 1 和 rho2\rho 2 是 MOD09A1 图像的蓝光和红光波段的图像反射率 (Yang et al., 2020)。重采样至 1-km1-\mathrm{km} 分辨率。使用最小-最大归一化将数据缩放到 0-1000-100 范围内。
SI=sqrt(rho_(1)xxrho_(2)), SSI=sqrt((NDVI-1)^(2)+SI^(2))
The rho1 and rho2 are the image reflectances of the blue and red light band of the MOD09A1 image (Yang et al., 2020). Resampling to 1-km resolution. Scale data into the 0-100 range using Min-Max Normalization| $\begin{aligned} & \mathrm{SI}=\sqrt{\rho_{1} \times \rho_{2}} \\ & \mathrm{SSI}=\sqrt{(\mathrm{NDVI}-1)^{2}+\mathrm{SI}^{2}} \end{aligned}$ |
| :--- |
| The $\rho 1$ and $\rho 2$ are the image reflectances of the blue and red light band of the MOD09A1 image (Yang et al., 2020). Resampling to $1-\mathrm{km}$ resolution. Scale data into the $0-100$ range using Min-Max Normalization |
Data set Year and format Data source Data processing Weight
Fractional vegetation coverage (FVC) " 2000, 2005, 2010,
2015, 2020; Raster, " 1km NDVI data: Juntao (2022) FVC =( NDVI - NDVIsoil )//( NDVIveg - NDVIsoil). Select the cumulative frequency of 5% and 95% of NDVI value as NDVIsoil and NDVIveg 0.16
Biological richness index (BRI) "2000, 2005, 2010,
2015, 2020; Raster,
1 km" Land use and cover data: www.resdc.cn BRI=Abioxx(0.35 xx forest +0.21 xx grassland +0.28 xx water area +0.11 xx cropland +0.04 xx construction land +0.01 xx unused land)/area. Abio is a normalized coefficient with a reference value of 511.264. Scale data into the 0-100 range using Min-Max Normalization 0.19
Net primary productivity (NPP) 2000, 2005, 2010, 2015, 2020; Raster, 500 m https://lpdaac.usgs. gov/products/ mod17a3hgfv061/ Resampling to 1-km resolution. Scale data into the 0-100 range using Min-Max Normalization 0.15
Soil salinization index (SSI) 2000, 2005, 2010, 2015, 2020; Raster, 500 m MOD09A1 image: https://lpdaac. usgs.gov/products/ mod09a1v006/ "SI=sqrt(rho_(1)xxrho_(2)), SSI=sqrt((NDVI-1)^(2)+SI^(2))
The rho1 and rho2 are the image reflectances of the blue and red light band of the MOD09A1 image (Yang et al., 2020). Resampling to 1-km resolution. Scale data into the 0-100 range using Min-Max Normalization" -0.14
Soil erosion 2000, 2005, 2010, 2015, 2018; Raster, 1 km https://www.resde.cn/ Scale data into the 0-100 range using Min-Max Normalization -0.13
Water conservation 2000, 2005, 2010, 2015, 2020; Raster, 1 km Wang (2022) Scale data into the 0-100 range using Min-Max Normalization 0.12
PM_(2.5) concentration 2000, 2005, 2010, 2015, 2020; Raster, 1 km https://sites.wustl.edu/ acag/datasets/surface-pm2-5/#V5.GL. 02 Scale data into the 0-100 range using Min-Max Normalization, and the maximum value is the average of the top 50 values in the 3 years -0.11| Data set | Year and format | Data source | Data processing | Weight |
| :---: | :---: | :---: | :---: | :---: |
| Fractional vegetation coverage (FVC) | $\begin{aligned} & \text { 2000, 2005, 2010, } \\ & \text { 2015, 2020; Raster, } \end{aligned}$ $1 \mathrm{~km}$ | NDVI data: Juntao (2022) | FVC $=($ NDVI - NDVIsoil $) /($ NDVIveg - NDVIsoil). Select the cumulative frequency of $5 \%$ and $95 \%$ of NDVI value as NDVIsoil and NDVIveg | 0.16 |
| Biological richness index (BRI) | 2000, 2005, 2010, <br> 2015, 2020; Raster, <br> 1 km | Land use and cover data: www.resdc.cn | $\mathrm{BRI}=\mathrm{Abio} \times(0.35 \times$ forest $+0.21 \times$ grassland $+0.28 \times$ water area $+0.11 \times$ cropland $+0.04 \times$ construction land $+0.01 \times$ unused land)/area. Abio is a normalized coefficient with a reference value of 511.264. Scale data into the 0-100 range using Min-Max Normalization | 0.19 |
| Net primary productivity (NPP) | 2000, 2005, 2010, 2015, 2020; Raster, 500 m | https://lpdaac.usgs. gov/products/ mod17a3hgfv061/ | Resampling to 1-km resolution. Scale data into the 0-100 range using Min-Max Normalization | 0.15 |
| Soil salinization index (SSI) | 2000, 2005, 2010, 2015, 2020; Raster, 500 m | MOD09A1 image: https://lpdaac. usgs.gov/products/ mod09a1v006/ | $\begin{aligned} & \mathrm{SI}=\sqrt{\rho_{1} \times \rho_{2}} \\ & \mathrm{SSI}=\sqrt{(\mathrm{NDVI}-1)^{2}+\mathrm{SI}^{2}} \end{aligned}$ <br> The $\rho 1$ and $\rho 2$ are the image reflectances of the blue and red light band of the MOD09A1 image (Yang et al., 2020). Resampling to $1-\mathrm{km}$ resolution. Scale data into the $0-100$ range using Min-Max Normalization | -0.14 |
| Soil erosion | 2000, 2005, 2010, 2015, 2018; Raster, 1 km | https://www.resde.cn/ | Scale data into the 0-100 range using Min-Max Normalization | -0.13 |
| Water conservation | 2000, 2005, 2010, 2015, 2020; Raster, 1 km | Wang (2022) | Scale data into the 0-100 range using Min-Max Normalization | 0.12 |
| $\mathrm{PM}_{2.5}$ concentration | 2000, 2005, 2010, 2015, 2020; Raster, 1 km | https://sites.wustl.edu/ acag/datasets/surface-pm2-5/#V5.GL. 02 | Scale data into the 0-100 range using Min-Max Normalization, and the maximum value is the average of the top 50 values in the 3 years | -0.11 |