陆地表层生态系统受到气候和人类活动的双重影响, 气候虽对生态环境变化存在显著响应, 但长时间尺度内趋于稳定, 难以发生大的突变^[1]。21世纪以来, 随着“人类世”概念的提出^[2], 人地关系日益紧密, 人类活动已然成为生态环境演变的重要驱动力之一^[3—4], 对自然环境的作用不容忽视。人类活动强度(Human Activity Intensity, HAI)是表示人类活动对自然环境开发利用程度的客观指标^[1]。人类足迹指数^[5—6]及影响力赋值法^[7—8]制备的HAI数据集借助空间模拟技术^[9—12]实现了HAI空间化表达, 但由于人类活动趋于复杂及多元化, HAI的空间化仍在发展, 且统一评价标准的缺失使各方法所得数据无法进行横向比较, 一定程度滞后了HAI数据的普适化。长时间序列与大空间尺度条件下量化人类活动面临较大的挑战, 但对深入探究复杂人地关系及其成因具有重要意义^[13]。

植被是生态环境变化的综合指示器^[14], NDVI时空演变已形成相对完善的研究体系, 多利用MODIS NDVI等数据^[15], 集中于黄土高原^[16—17]、青藏高原^[18]等生态保护区。结果显示, NDVI对气候因子(降水、气温等)和人类活动存在响应^[19], 在生态脆弱区反应更为敏感, 但对人类活动响应机制缺乏深入认识。人类活动强度研究也多聚焦单因子时空分异^{[16, 20]}, 与自然因子联动分析较少, 多借助真实值与拟合值的残差^[17]表示人类活动对环境的扰动。事实上, 人类活动对自然环境的作用具有双向性^{[16, 21—22]}：城市扩张、资源开采等活动会明显扰动原有生态系统运转, 表现为植被生物量减少, 红光波段反射性增强, NDVI减少, 二者显著负相关^[23]；自然保护区建设、封山育林、开荒造林等生态恢复策略对区域生态环境改善与维持稳定有显著正向驱动^[24]。人地相互作用是多要素共同决定的复杂过程, 植被退化与修复取决于人类活动的方向性, 且存在强烈区域差异及阶段性, 人类活动对生态环境综合影响的研究亟待完善。且由于人类活动高频区域的空间分布具有局限性, 大尺度下其影响易被忽略, 中、小尺度下人类活动对植被更具影响力^[16]。目前中、小尺度相关分析犹显不足, 选取合适尺度是探究二者时空关联的关键。

黄土高原位于西北干旱区平原向高原的过渡地带, 坡面物质不稳定, 水土流失严重, 是国家生态修复工程重点实施地。区域内生境质量、人类活动的相关研究多从植被覆盖时空演变、土地利用变化和人类活动强度等单方面展开, 涉及区域HAI与NDVI的关联性研究还有待深入。故基于土地利用、夜间灯光及人口密度数据构建人类活动强度表征模型, 采用斜率分析法对黄土高原人类活动强度与植被覆盖时空分异进行分析, 并用相关分析、空间自相关及标准差椭圆等方法探究二者的时空关联, 旨在揭示区域人地关系发展态势, 以期为控制资源开发利用程度、调节人类活动强度、维护生态环境与风险防范提供参考。

1 研究地区与数据来源 1.1 研究区概况

黄土高原是四大高原之一, 位于中部偏北地区(33°39′—41°16′N, 100°52′—114°31′E)^[25], 地势西高东低(图 1), 主要有山西高原、陕甘晋高原、陇中高原、鄂尔多斯高原和河套平原组成, 面积约为62.60万km^2[26], 以其黄土覆盖面积广及地势较高得名。为半湿润半干旱过渡区, 气温及降水自东南向西北呈带状分布且逐级递减, 年平均降水量在200—700mm^[15]。黄土高原地跨山西省和宁夏回族自治区全域、陕西省北部、内蒙古自治区、青海省和河南省部分地区在内的7省45地级市, 形成了一个区域性城市群(关中平原城市群)和四个地区性城市群(呼包鄂榆城市群、晋中城市群、宁夏沿黄城市群和兰西城市群), 经济发展已形成较成熟格局^[27]。由于黄土性质特殊, 汛期多暴雨, 人类活动加剧, 水土流失问题突出, 人地矛盾尖锐^[28—29], 21世纪后进入植被大规模、快速恢复和重建时期^[30—31], 人类活动强度与植被覆盖间存在一定时空关联性。

1.2 数据来源

研究采用DMSP/OLS(2000—2013年)、NPP/VIRS(2012—2020年)两类夜间灯光数据, 来源于美国国家海洋和大气管理局的国家地球物理数据中心。由于VIIRS数据传感器过境时间相对延迟, 需将两种数据源统一投影, 并分别校正和整合校正^[32], 对DMSP/OLS数据进行饱和和连续性校正, 对NPP-VIIRS数据进行背景噪声与异常值消除^[33]。以2013年校正后的DMSP/OLS数据为基准, 对两种数据进行拟合, 选取最优回归模型进行2014—2020年NPP-VIIRS向DMSP/OLS的转换, 得到类DMSP/OLS数据集。DEM和土地利用数据由中国科学院资源环境科学与数据中心(https://www.resdc.cn/)获得, 人口密度数据来源于WorldPop网站(https://www.worldpop.org/), 2000—2020年生长季(5—9月)归一化植被指数来源于国家青藏高原科学数据中心(https://data.tpdc.ac.cn/zh-hans), 并利用最大值合成法^[34]获得2000—2020年逐年NDVI数据。所有数据分辨率为1km, 投影为WGS_1984_UTM_Zone_48N。

2 研究方法 2.1 人类活动强度表征

人类活动多元化、复杂化致使统一量化标准难以确定, 前期应用的是基于统计数据的非空间化方法^{[9, 35—36]}, 后期空间化概念开始形成^[15]。其中, 以土地利用类型变化法^[37—39]、综合指标法^{[12, 40]}及人类足迹指数^[6]为代表。空间化方法虽选取土地利用、人口密度、夜间灯光数据^[8]等计算影响力进行评价, 但精确度较差, 无法细化区域内部差异。结合徐勇^[37]、陈泓瑾^[41]等人的研究成果, 以土地利用数据为基础, 综合陆地表层人类活动强度(HAILS)、夜间灯光及人口密度数据, 构建HAI表征模型：

式中, HAI人类活动强度, NTL、PD、HAILS为归一化后的夜间灯光、人口密度、陆地表层人类活动强度数据, S_CLE为建设用地当量面积, S为区域总面积, SL_i为第i种土地类型面积, CL_i第i种土地类型建设用地折算系数, n为区域内土地类型数目。参考相关研究^[41], a、b、c分别取值为0.3、0.3、0.4, 建设用地、耕地、草地折算系数^[20]分别为1、0.4、0.067, 其余地类折算系数为0。

2.2 线性回归法

利用一元线性回归法进行NDVI及HAI长时间序列变化趋势分析^[42], 公式如下：

式中, Slope为NDVI与HAI变化速率, n为年数, 本文取21；i为时间序列, 2000—2020年依次取1—21；α_i为各年NDVI及HAI值；Slope值正负表示NDVI/HAI随时间增/减, 其绝对值表示变化速率。

2.3 相关系数

相关系数用以表示两个变量线性相关程度, 取值区间为[-1, 1], 公式如下：

式中, r为NDVI与HAI相关系数；n为时间序列, 取值1—21；a_i、b_i分别表示归一化的NDVI与HAI值, r值正负表示二者呈正相关或负相关, 绝对值表示相关程度。

2.4 双变量空间自相关

使用双变量局部空间自相关分析黄土高原NDVI与HAI空间关联关系、空间聚集区及其分布^[43], 公式如下：

式中, W_ij为空间单元i与j的空间权重邻接矩阵(i, j=1, 2, 3, …, n), x_i、y_i为观测值；为x_i的平均值；I取值[-1, 1], 利用p值进行显著性检验。

2.5 标准差椭圆及重心模型

标准差椭圆和重心计算可获取高植被覆盖与高人类活动强度区域重心及其时空移动轨迹, 可体现其空间形态演变、主要驱动方向及要素空间分布的集散性^[44]。

2.6 HAI变化频率图谱构建

参考构建土地用变化图谱相关方法^[45], 将栅格尺度下黄土高原HAI空间分布及其演变过程借助编码重现, 结合图谱代码融合方法^[46], 选取六个时段, 获得六位编码的HAI时空演变图谱和长时间序列HAI变化频率, 频率图谱栅格属性值取值为0—5。

3 结果分析 3.1 黄土高原人类活动强度时空特征

黄土高原2000—2020年HAI年均变化较缓, 变化率为1.39%/a。空间差异显著, 根据自然断点法分为低、较低、中、中高、高强度五个等级, 人类活动强度整体空间分布呈东南高西北低(图 2)：中高、高强度区集中位于河谷平原区等低海拔地带, 占比6.29%；中强度区占比17.06%, 位于中高强度区外围；低及中低强度区占比76.65%, 位于中部丘陵、高塬沟壑区及西北沙地和农灌区。

研究时段内HAI整体呈现“七点一圈一带”圈层分布模式(图 2), “七点”为黄土高原所辖省会和发展程度高的城市, 是人类活动的中、高强度分布区, 逆时针排序依次为：西安市、太原市、大同市、呼和浩特市、包头市、银川市及兰州市, 城区HAI达到0.70以上；“一圈”为七点城市联通构成外围中高人类活动强度圈层, 途径关中城市群、晋中城市群、呼包鄂城市群、宁夏沿黄城市群及兰西城市群, 在省会城市的辐射作用下, 人类活动强度达到中高等级, 以中心建成区为高强度区向外辐射, HAI变化遵循距离递减规律^[33]；中部呼和浩特市以南, 鄂尔多斯市以东, 吕梁山脉以西, 天水市以北, 自东北向西南部形成狭长中低强度人类活动带。2017年后带状形态更为显著, 依次经过榆林市、延安市、庆阳市、平凉市、固原市、天水市等, 受地形地貌与距核心城市较远影响, 辐射作用较弱, 平均HAI为0.21。

低HAI区域分布位置较为特殊, 多为生态保护区及沙漠等未利用地, 具体位于延安市子午岭自然保护区、晋西太行山脉、南部秦岭山脉, 陇、宁及鄂尔多斯市内的沙漠地区。以西安市为代表, 高强度城市HAI变化率呈现圈层式变化, 中心城区-城市扩张区-外围乡镇HAI变化率为低-高-低。原有建成区活动强度变化不大, 增长率不超过0.27%, 扩张区域HAI增长速率较快, 变化率在1.04%—4.60%, 周边农村居民点人类活动强度以相对较低速率(0.27%—1.04%)增长。中强度以上区域HAI变化较显著, 低强度区域受地理条件限制, 2000—2020年变化速率较低, 难以进行多样人类活动。

长时间尺度下, 黄土高原人类活动强度变化频率空间分异明显。以2020年HAI值为基准, 4年为节点, 划分人类活动强度, 获得基于栅格尺度的2000—2020年黄土高原HAI变化频率图谱(图 3), 将人类活动强度变化频率划分为稳定不变区(HAI转换次数为0)、低频转换区(转换次数为1—2)、高频转换区(转换次数为3—5)。其中, 稳定不变的区域占比71.71%, 用地类型以林地、未利用地(沙地)居多, 活动强度小, 土地利用结构稳定。低频转换区覆盖面积总占比24.06%, 集中分布于原有建成区的外围地带, 已有一定人类活动基础, 长时间序列下不会发生过高频率的转化, 表现为城市建成区为中心向外扩散。以西安市、太原市为代表, 人口逐渐向外缘地区聚集, HAI呈现圈层式增长, 同时反映出2000—2020年黄土高原城市呈现扩张态势, 与已有研究结果一致^[39], 在河套平原等区域人类活动也有小幅变化。

高频变化区面积总占比不足5%, 在鄂尔多斯市与榆林市交界处等地聚集性分布(图 3), 与高强度人类活动区域分布位置吻合(图 2)。该区域地处呼包鄂榆城市群中心位置, 是内蒙古自治区和陕西省经济发展的重要增长极^[47]。自确立以呼包鄂为核心的特色经济圈发展战略和将呼包鄂榆城市群纳入全国城市群规划建设后^[48—49], 该区域经历了明显的城市扩张, 形成了呼市为核心, 包鄂榆三市为增长极的发展模式, 且发展重心向榆林市移动。城市联动明显增强, 以鄂尔多斯市东胜区和榆林市榆阳区最为典型, 人类活动强度实现高频转化, 体现了人类活动与城市发展的同向性。

3.2 黄土高原NDVI时空特征

2000—2020年黄土高原植被NDVI呈现波动上升趋势, 年均增速为1.59%, 生态恢复效果显著。2000— 2013年增速较快, 实现每年2.11%的增速上涨, 2014—2020年增速有所放缓, 以每年0.82%的速度提升。生态工程对植被恢复的贡献率可能符合环境库兹涅茨倒U型曲线, 即在一定时间内退耕还林还草对改善生态环境起到正向作用, 但作用效果先升后降, 达到最大值后, 其影响将有所减弱。后期在推进生态建设的条件下, NDVI增速减缓可能与植被自然演替过程有关。

受气候因子及人类活动等作用^[50—51], 黄土高原植被NDVI及其变化率存在明显空间分异：以自然断点法将NDVI年均值分为低、较低、中、中高及高植被覆盖等级, 自东向西依次为中高覆盖带、中低植被覆盖带、低植被覆盖带, 阶梯型分布显著(图 4), 且中高植被覆盖区域向西北扩张(面积占比42.48%), 集中体现于陕北及甘肃省部分地区, 低覆盖区域显著缩小(面积占比14.63%)。NDVI变化率自东南向西北呈低-高-低形态：关中城市群与晋中城市群发展程度高, 建设用地涨幅大, NDVI增长速率较慢, 西安、咸阳等市出现负增长, 达到负向变化率最高值3.57%；中部陕北以及陇东部分地市是国家黄河中上游重点生态环境建设区^[52—53], 同时作为退耕还林还草重点实施区^[54], 在东北西南方向形成NDVI高速变化带, 实现最大增长率4.45%, 植被恢复效果显著；鄂尔多斯市、宁夏回族自治区以及甘肃省部分地市受土地类型及地势地貌影响, 植被恢复难度高, 速度慢。总体来看, 植被恢复区域占比高达96.32%, 植被覆盖负增长区域占比低于4%, 黄土高原整体植被情况向好。

3.3 黄土高原NDVI与HAI时空关联性分析

由于黄土高原人类活动区域相对固定, 结合图 4分析可知NDVI与HAI时空关联性较强区域与主要人类活动区相吻合(图 5)。根据相关系数大小分为显著正相关、较显著正相关、不显著相关、较显著负相关、显著负相关五类。HAI与NDVI均降低使二者呈正相关的区域占比不足0.5%, 故本文只讨论HAI增加的同时, 植被覆盖得到明显缓解的正相关区域。该区域面积占比达19.39%, 显著正相关区占比1.35%, 集中分布于呼包鄂榆城市群及延安市、庆阳市内, 其余地市零星分布, 负相关区域占比9.54%, 显著负相关区域面积覆盖比为1.42%。由于人类活动有双向作用, HAI与NDVI呈现负相关有两种形成机制：一是资源消耗型负相关, 人类活动强度增加, 用地类型以建设用地和耕地为主, 植被覆盖率较低；二是资源补偿型负相关, 人类活动强度下降, 退耕还林还草以提高植被覆盖。以2000—2020年NDVI变化率为负区域(图 5)为基准, 将负相关区域分为上述两类, 消耗型负相关区域基本与城市建成区位置吻合, 城市建设导致人类活动加剧, 植被覆盖情况恶化。城市周边由于建成区外扩及农村居民点规模增大, HAI加强, 负相关显著性较低, 仍属于消耗型负相关区域。补偿型负相关区域集中位于山西省朔州市、陕北榆林市、延安市、以及甘肃省平凉市, 为NDVI高速变化带(图 4), 在其余区域零星分布。

3.3.1 双变量空间自相关

为进一步探究NDVI与HAI时空关联性, 依托双变量空间自相关模型, 以10km格网为单位, 得到2000—2020年二者空间关联模式(图 6)。H-H表示高强度与高覆盖聚集；L-L表示低强度与低覆盖聚集；L-H表示低强度与高覆盖聚集；H-L表示高强度与低覆盖聚集, N为分布无显著聚集特征区域。空间上以不显著区域为界, 以东为高强度-高覆盖和低强度-高覆盖区域, 以西为低强度-低覆盖和高强度-低覆盖区域, 分异明显。

高高聚集区分布于东部关中及晋中城市群, 低高聚集区与高高聚集区相邻分布, 与秦岭山脉、子午岭及吕梁山脉位置基本吻合, 用地类型多为林草地, 植被覆盖情况良好；低低聚集区在西北方向片状分布, 内辖库布齐沙漠和毛乌素沙地以及腾格里沙漠的部分区域, 植被生长条件差, 人类活动强度低；自南向北沿黄河流域部分区域人类活动强度较高, 高低聚集区在低低聚集区外围零星状分布。

时间尺度上, 高高聚集区、低低聚集区、高低聚集区均呈现收缩态势, 面积分别下降0.24%、3.79%、1.90%, 高强度区对生态环境的扰动得到了有效控制, 不适宜植被生长的区域生态恢复效果显著, 高低聚集区人类活动强度与植被覆盖差距有所缓解；低高聚集区呈扩张趋势, 主要由高高聚集区及不显著区转换而来(图 6), 原有林草地得到了妥善保护, 退耕还林等修复工程取得了鲜明成效, 高强度人类活动区蔓延态势得到有效抑制。不显著区变化幅度最大, 面积同比增长4.30%, 相邻区域差异缩小, 表明黄土高原整体人类活动与生态环境向协调方向发展。

3.3.2 重心转移与标准差椭圆

借助重心模型与标准差椭圆描述各期高植被覆盖区(NDVI≥0.6)与高人类活动强度区(HAI≥0.6)时空移动轨迹(图 7)。2000—2020年黄土高原NDVI重心始终在陕北榆林与延安西部交界处, 2000—2010年向西南移动, 2010—2020年向东北移动, 整体向西北转移。2000年开启退耕还林等工程以来, 延安市及陇东城市群是重点实施区域^[52], 整体植被恢复效果优于北部鄂榆城市群的沙漠地带, 后期鄂榆地区库布齐沙漠、毛乌素沙地等区域的生态工程对NDVI重心产生了一定拉力效应。标准差椭圆面积自2000年39.43万km²增长至2020年43.17万km², 表明高植被覆盖区域向西北方向增大且多方向外延, 修复效果显著。标准差椭圆扁率无显著变化, 多年均值为24.08, 高覆盖区总体仍沿东北-西南方向倾斜分布。

2000—2020年高强度活动区重心显著北移, 在东西方向上稍有转移, 整体向西北移动, 西北区域HAI增加, 发展水平显著提升。2000—2015年重心位置相对集中, 年均移动速度为2.50km/a, 2015—2020年年均移动距离为12.60km, 北移距离显著高于西移。已有研究表明, 2010—2018年间, 呼包鄂榆城市群大规模扩张, 增长率为97.70%^[53]。依照规划将建设成为中西部地区核心城市群^[49], 同时宁夏沿黄城市群形态更为成熟。西北城市的迅速崛起, 拉动该区人类活动强度大幅提升, HAI高值区重心向西北方向移动。标准差椭圆面积由2000年44.48万km²减少至2020年38.09km², 同比降低14.36%, 表明高强度人类活动区域逐渐集中。椭圆扁率在北偏东方向先增后减, 前期(2000—2010年)东南区域发展较快, 高强度区朝东南方向集中, 后期(2010—2020年)西部城市群发展加快, 开始向西北转移。高覆盖与高强度区重心均朝西北方移动, 侧面解释前文黄土高原中部区域HAI与NDVI呈现正相关, 同时表明生态保护与人类活动空间上显现出同向性, 人地关系协调程度提高。

4 讨论

整体看来, 黄土高原人类活动与植被覆盖时空关联性呈现“北正南负、西正东负”的空间格局, 负相关区域主要集中在关中平原城市群、晋中城市群和陕北两市交界处；与植被恢复工作协调性整体呈现“北高南低、西高东低”的时空特征, 重点恢复区为陇东黄土高原、宁夏平原及西北防护林带(图 5、图 8)。以延安市为例, 分析正相关区域的背后成因。延安市地处黄土高原腹地, 分布有大量煤炭、天然气等矿产资源^[55], 是典型的资源导向型城市, 位于黄土丘陵沟壑区之上, 也是我国代表性生态脆弱城市^[56]。1999年和2013年先后开展了两轮大规模退耕还林还草、封山育林等工作, 保水保土成果丰硕, 中北部植被覆盖大幅好转(图 8), 区县尺度上以北部子长县、安塞区及西部吴起县、志丹县效果最为明显^[55]。北部植被恢复速度整体高于南部, 土地类型以耕地向林地、草地的转入居多, 侧面佐证了前文黄土高原高植被覆盖区域重心北移的迹象。人均定量基本农田建设、大力发展经济林、城市化和生态移民等系列政策^[57—58]的相继出台增加了农民收入, 城镇人口增长迅速, 西部大开发后社会经济快速发展^[59], 人类活动强度显著提高。研究表明, 近20年延安市土地利用结构趋于稳定, 优势土地类型主导作用增强^[60]。2000—2010年, 北部子长、吴起、志丹三县受经济与生态工程影响, 地类转换频繁；2010年后生态治理由规模扩张向成果巩固阶段转化, 经济发展也由高速增长向高质量转化^[61], 土地利用结构趋于稳定, 建设用地的增长及竞争优势均高于耕地, 且呈扩张趋势, 解释了延安市西北区域2000—2020年人类活动强度正向变化的形成。其他正相关区域也受到人类活动的双向影响：呼包鄂榆城市群位于包昆经济带沿途, 是我国“两横三纵”经济带要素畅通的重要载体^[62], 同时毛乌素沙地和陕北是我国三北防护林区域^[63]和植树造林、退耕还林重点实施区, 人类活动强度与植被覆盖同步提高；晋北朔州、大同等地由于开采矿产^[64]、城市扩张等原因建设用地激增, 得益于天然林保护和黄土丘陵沟壑区整理工程实施, 植被覆盖显著好转^[19]。

黄土高原位于国家“两屏三带”生态安全战略区, 同时地处黄河中上游, 是黄河流域生态保护和高质量发展国家战略的关键实施区域, 人类活动要兼顾经济与生态成效。东南区域发展程度高, 建设用地与耕地广布, 西北分布有三北防护林、防风固沙区和水土保持区, HAI的空间差异性和生态建设工程作用下黄土高原形成了“西密东疏”的生态网络结构^[65]。资源消耗活动虽对植被生长与生态保护造成一定扰动, 但开展补偿修复工程可一定程度缓解环境恶化。为保证区域生态环境和高质量发展, 应结合土地利用结构制定差异化政策, 东部在稳定经济发展基础上加大力度保护现有生态源地, 建立绿化修复区, 建成区内实行小范围生态工程, 同时优化西部屏障区的生态结构, 把握好生态环境建设与社会经济发展的平衡。

5 结论

综合夜间灯光、人口密度、土地利用数据量化表征人类活动强度, 基于栅格尺度分析2000—2020年黄土高原HAI、NDVI时空演化特征, 二者的时空关联性及其空间分异, 结论如下：

(1) 20年间黄土高原HAI时间尺度变化较平缓, 空间尺度分异较大, 呈“七点一圈一带”的分布模式, 呼包鄂榆城市群及晋中城市群人类活动强度等级变化频率高。NDVI整体增速为1.59%/a, 空间尺度始终呈东南高西北低, 陕北等退耕还林实施区植被覆盖度显著改善。

(2) 2000—2020年黄土高原高强度人类活动与高植被覆盖区重心均向西北方偏移, 高强度人类活动区蔓延态势得到有效抑制, 植被贫瘠区生态恢复效果显著, 人类活动与植被恢复空间呈现同向性, 人地矛盾缓解, 协调发展程度提高。

(3) 人类活动对植被覆盖有双向作用：时空关联性呈现“北正南负、西正东负”的空间格局, 东南经济发达区呈负相关, 西北生态修复区植被覆盖明显改善, 区域整体对生态恢复的积极效应逐渐提升, 能够较好把握生态保护与高质量发展的平衡。