自进入“人类世”以来, 人口的迅速增长以及人类对地球资源的开发利用, 对全球自然环境及生态系统造成了严重影响^[1]。由于特定自然环境和人文条件对人口全空间布局及人类有序活动的限制^[2], 全球范围内仍有大量未被人类改造、或只被轻微改造但保留有效生态功能的“荒野地”^[3—4], 这些区域能够提供重要的生态系统服务功能^[5—6], 在生物多样性保护^[7]、抵御气候变化^[8]、改善人类福祉实现可持续发展^[9]等方面具有重要作用。目前, 许多国家和地区已经通过法律或行政手段建立了荒野地保护机制, 并且进行了大量的实践^[10—11]。中国尚未在法律和政策中直接提及荒野地概念^[12], 学术界关于中国荒野地的研究也才刚刚起步, 面向中国陆域开展荒野地研究工作, 能够为完善中国现行生态环境治理体系提供参考, 对自然生态空间保护工作具有重要意义。

荒野地空间识别是荒野地研究的首要任务, 能有效揭示荒野地的空间分布格局、挖掘其时空演变规律^{[2, 13]}。现有的荒野地识别模型多基于多源空间数据建立指标体系, 利用多指标评价方法(Multi Criteria Evaluation, MCE)进行指标叠加, 实现对荒野地空间范围的识别^{[10, 14—15]}。大尺度荒野地研究面向国家^[12]、大洲^[16]乃至全球^[17], 注重对荒野地分布格局及演变规律的探究^{[2, 18]}, 其指标体系的构建依赖于数据源的可获得性。例如, 荒野连续谱法强调自然环境由城市端向自然端的过渡, 认为一个地区更加远离人类干扰, 则其更可能是荒野地^[15], 通常使用自然度、遥远度、崎岖度、破碎度等表征自然水平和偏远程度的空间指标^[19]; 人类足迹法基于人口密度、道路、铁路、农田、牧场、建筑区和夜间灯光等数据, 建立人类压力指标以量化人类对生态环境的影响^[17—18]。而高精度时空数据支持下的小尺度荒野地研究则聚焦区域^[20], 以及单个保护地^[19], 更加关注所提取边界的准确度, 评估面向复杂空间特征的荒野程度, 其指标体系的构建更加精细且具有地域特点。例如Brackhane等^[11]考虑生态、地理和法律政策等多个方面, 基于生境连续性、外来物种、噪声污染等二十个指标识别了Rechlin德国国家自然遗产公园的荒野地范围。在指标叠加方面, 目前常用的方法是线性加权组合(Weighted Linear Combination, WLC), 对多维度多类型指标归一化后加权叠加, 得到一个可表征荒野程度的连续变量, 然后在荒野度区间内确定符合荒野地条件的阈值^{[10, 15]}, 实现对荒野地的识别。

以往的研究多采用等权重相加的指标组合方式, 由于指标之间的相关关系和重要程度难以确定, 可能会低估或者高估某个指标的影响。另外, 荒野地本质上是人类构建的主观概念^[21], 荒野和非荒野之间分界点的确定并没有统一的科学标准, 可能导致识别结果存在一定程度上的差异^[15]。荒野虽然是个相对概念, 但其对立面具有明显界限, 例如城市、农田、牧场、现代交通路网等是被广泛认知的非荒野区, 在明确非荒野地的基础上量化潜在荒野地的荒野程度, 进一步识别荒野地的空间范围, 可为荒野地评估提供具有实践价值的视角和方法。例如, 马力等^[2]建立了一种决定—影响的耦合关系模型, 通过气候和地形因子确定荒野地分布范围, 选取土地利用、夜间灯光、居民点和道路因子确定荒野地质量。

综合以上原因, 本文提出一种适用于中国陆域荒野地识别的方法模型, 构建生物物理自然度、可达性、人口密度、人工环境和光污染影响五个指标; 纳入高分辨率空间数据, 考虑指标对象的地理特征, 通过构建与精度相适应的降尺度算法和指标计算模型提高了识别精度; 针对传统WLC方法的局限性, 考虑面向复杂地理环境的人类活动模式, 构建多要素综合的计算模型, 并利用加权乘积和布尔运算进行指标叠加分析, 提供了一种可实践的指标综合思路。

1 数据来源与研究方法 1.1 数据来源

本文用作计算的数据包括DEM、土地覆盖、居民区、交通路网、人口密度、牲畜分布和光污染数据。将所有栅格数据重采样为30m×30m, 为获取较高精度的评估结果, 进一步基于DEM、土地覆盖数据等对牲畜分布数据作降尺度计算, 并利用3×3的窗口对人口密度和光污染数据做均值滤波, 以平滑低分辨率向高分辨率转换过程中的数值突变。基于人类足迹法的荒野地分布数据用作与本文的荒野地识别结果进行对比实验。数据来源及相关信息见表 1。

1.2 整体思路

本文参考世界自然保护联盟对荒野地的定义^[15], 将其描述为无明显人类影响或只被轻微影响、仍保留其自然特征的区域, 然后利用现有数据和GIS分析方法, 定量研究人类活动对生态的影响强度并构建荒野地识别指标体系。考虑源数据的可计算性, 将人类对自然的影响主要概括为: 人类非持续性活动对自然的干扰(通勤、劳作、徒步旅行、狩猎等人类出行行为)、人类长期改造自然形成的人工环境(农田、牧场、建筑区和道路)以及夜间人造光对生物环境的光污染影响^{[2, 18, 28]}。基于此, 构建荒野地识别指标体系及计算方法: (1)分别建立生物物理自然度指标反映面向不同土地利用类型时人类出行选择的偏向性、可达性指标表征人类出行的空间范围、人口密度指标量化出行的数量特征, 基于此三个指标利用等权重相乘模拟人类非持续性活动强度; (2)分析人类长期改造背景下的地表生态景观演化结果, 建立人工环境指标, 以模拟人类长期或永久性改造活动的空间分布格局; (3)引入新型光污染数据建立光污染影响指标, 以识别超越人类夜间活动物理边界的光污染影响范围; (4)基于人工环境和光污染影响指标确定人类长期活动区, 利用空间布尔叠加识别潜在的荒野地, 并综合人类非持续性活动强度, 得到国家尺度的荒野度分布数据, 具体流程如图 1。

1.3 指标计算 1.3.1 生物物理自然度

生物物理自然度反映了人类定居、农业活动等行为对生态系统原始状态的改造程度, 描述了不同土地利用类型的自然水平, 一般通过专家打分的方式来定量评价。本文基于中国土地利用数据(1980—2015)的分类体系^[29]对土地覆盖类型进行类型分类和综合, 并参照Cao等^[15]的标准评价生物物理自然度(表 2)。

1.3.2 可达性

可达性能有效表征人类出行范围, 是定量评价人类活动强度的可靠因子^[30—31]。过去的研究常基于居民点或道路网络, 使用欧式距离衰减法来量化可达性的分布范围, 无法全面反映复杂地理环境对实际人类可达性分布的影响。本文将可达性定义为, 某空间区域邻近居民到达此区域的最小累积时间耗费。具体来说, 考虑人类出行对多类型地表覆盖、交通路网和地形的响应特征, 引入成本距离算法度量由居民区到达邻近区域的最小累积时间耗费, 构建区别于以往欧式距离量化人类活动范围的计算模型, 以表征人类出行活动的空间特征。

地表覆盖类型和交通路网作为地表景观类型, 与出行特征具有高耦合性, 人类在不同地表景观内的移动速度是影响出行规律的重要因素。因此, 构建基于地表景观类型的通行速度指标, 并建立地表景观类型与通行速度的映射关系: (1)基于Weiss等^[32]研究成果, 综合GLC_FCS30数据集的土地覆盖类型, 并赋予步行速度; (2)考虑高速公路和铁路一般只提供居民区到居民区的通道, 划分其为封闭交通线路并赋值为0;(3)除高速公路以外的开放式公路, 提供了人类向非居民区快速扩散的通道, 以中国一级公路最低限速60 km/h为其通行速度; (4)步道是步行的首选, 赋予人类最高的步行速度5 km/h^[32]。结果如表 3。

在宏观视角下, 地表景观更多的是在水平面上影响人类活动, 而坡度作为一种地形因子, 则是在垂直面上影响通行速度, 因此建立坡度影响因子指标。仅考虑坡度为正的情况时, 通行速度与坡度呈负相关^[32], Tobler等^[33]提供了坡度-速率函数(公式(1)), 以此计算坡度影响因子(公式(2)):

(1)

(2)

式中, F_s为坡度因子, V_slope即坡度为θ(度数)时的步行速度(km/h)。

成本距离算法(Cost Distance)一种Dijkstra最短路径算法的变体, 通过计算在成本栅格上移动的最小累积耗费, 度量由源区域到达其它区域的最小成本耗费。本文综合通行速度指标和坡度因子指标, 利用公式(3)建立通行阻抗成本栅格, 像元值代表在此单元格移动的时间耗费。选取居民区数据为源数据, 利用成本距离工具, 计算从居民区出发到达邻近非居民区的最小累积时间耗费, 以量化人类到达此区域的难易程度, 表征出行的空间特征, 得到可达性指标(图 2), 可达性数值越高, 代表到达此区域的时间成本越高, 荒野度越低。

(3)

式中, T为通行时间成本, 单位h; L为通过每个像元的距离, 因本文采用30 m分辨率栅格, 故取值0.03 km; V为通行速度, 单位km/h; F_s为坡度因子。

1.3.3 人口密度

人口密度是影响人类活动强度的重要因素, 人口数量的增加, 加速了自然资源的枯竭和生态环境的退化^[34]。人口密度被广泛用于表征人类活动强度, 也是荒野地识别的重要数据来源^{[18, 27]}。基于此, 本文建立人口密度指标, 为获得全范围的空间数据, 利用极小值0.0001填补NoData和值为0的区域。

1.3.4 人工环境

人类在长期的改造自然活动中, 建立了城镇、道路、铁路、农田和牧场等人工设施, 深刻改变了当地的原始自然景观并在其周围一定范围内可能造成持续性的影响, 此区域人类干扰程度高, 为典型的非荒野地区。基于土地覆盖数据的不透水面提取城镇和建筑区空间数据, 同时将旱作农田、灌溉农田和乔木或灌木覆盖(果园)三种类型的土地覆盖作为农田数据。

道路、铁路等现代交通路网的建造改变了其周边的物理化学条件^[35—36], 产生了超越自身空间范围的边缘效应^[37], 从而形成了包含交通路网及其外延区的人工环境。交通路网对自然环境的实际影响程度和范围取决于路网类型、当地环境^[38]和交通流量等^[35—36], 因此引入缓冲区分析工具建立道路、铁路影响区模型, 具体为: (1)公路及周边200 m区域内, 重金属等化学污染含量较高^[37], 且大型哺乳动物密度较低^[39—40], 800 m范围内部分鸟类、两栖类、爬行类等受到明显影响, 故将800 m作为面向道路的人工环境范围阈值; (2)铁路及周边600 m范围内, 栖息地物种数量明显减少, 为生态风险较高区域^[38], 铁路的人工环境范围阈值设为600 m; (3)与公路和铁路相比, 步道的基建水平不高, 且大型交通工具出现频率较低, 因此对生态环境产生的化学污染以及噪声、电磁等物理属性的改变较少, 本文将100 m作为面向步道的人工环境范围阈值。

牧场是人类干扰生态环境的直接因子^[41], 已被用来表征人类活动影响^[42], 也是荒野地识别的重要指标^[27]。Ramankutty等^[43]利用农业普查数据和土地覆盖数据建立了10 km分辨率的全球牧场分布数据集, 在荒野地研究领域被广泛应用^{[18, 27]}。本文引入精度相对较高的牲畜分布数据集, 其包含多种牲畜的空间分布和数量信息, 为模拟较高精度的草原牧场空间分布, 构建了一种数据降尺度方法, 具体思路: (1)从GLW2数据集选取牛、山羊和绵羊三种牲畜的分布数据, 按照5、0.8和1的权重^[44]进行栅格叠加计算, 生成新的牲畜分布栅格; (2)在青藏高原地区采用QTP_DA_2020_Livestock_Seasonal_Data数据集, 其基于牧场位置数据和调查数据^[24], 提供了精度更高的青藏高原冷暖季牲畜空间分布和数量信息; (3)考虑草原牧场的人文条件、地表类型和地形特征, 建立草原牧场的条件因子及分类标准(表 4), 基于以上牲畜分布栅格及土地覆盖、DEM和坡度数据进行逐栅格条件计算, 构建符合草原牧场条件的新栅格, 识别中国草原牧场的空间范围(图 3)。

1.3.5 光污染影响

光污染是人造光源引起的夜间自然光照水平的改变^[46], 可对生物过程、生物多样性和生态系统的功能产生负面影响^[47—48], 是人类对生态的重要干扰因子^[49—50]。人工夜空亮度数据基于VIIRS/DNB遥感影像和天空亮度测量数据计算, 能够定量描述光污染程度, 依据Falchi等^[26]研究成果, 人工夜空亮度与自然夜空亮度(174 μcd/m²)的比值可评估光污染的影响水平, 以此建立光污染影响指标(图 4)。其中, 指标大于0.5时, 表明环境已经受到明显的光污染影响, 天空的自然外观已经消失, 为人类干扰强度较高区域^[26]。

1.4 荒野度计算

人类非持续性活动主要反映为出行行为, 一般而言, 人类更倾向于选择容易到达和有吸引力的地方作为目的地。其影响强度主要取决于出行选择的偏向性、到达的难易程度和区域人口的数量特征, 分别以生物物理自然度、可达性和人口密度三个指标来表征此三个要素, 计算人类非持续性活动对生态的影响强度, 表达式如下:

(4)

式中, HAI为人类非持续性活动强度, 数值越大表明人类对自然的影响越大, 荒野程度越低; BN为生物物理自然度指标; HA为可达性指标, HA_max为可达性的数值最大值; PD为人口密度指标。

人工环境和光污染覆盖区域即人类长期活动区, 是人类影响程度最高的区域, 为典型的非荒野地区, 利用布尔运算将其与人类非持续性活动强度叠加, 然后按照人类非持续性活动强度越高荒野度越低的规则评价荒野度, 得到全国范围内的荒野度分布数据(图 5)。

2 结果与分析 2.1 荒野度分布格局

中国陆域荒野度分布具有典型的空间分异性, 呈现出“西高东低”的特征, 荒野度面积占比如表 5。荒野度较高的区域(3—5)主要分布在中国西部、北部(图 5), 荒野度较低的区域(1—2)主要分布在中国东部、中部和南部, 这些地区人口密度较高, 人类活动相对频繁。人类长期活动区为占全国面积的三分之一以上(37.56%), 主要分布于四川盆地、黄土高原和东北、华北、长江中下游等平原区域, 这些地区地势平坦, 适于人类居住, 且社会经济发达、交通便利, 人类活动最为频繁, 是典型的非荒野地。

2.2 荒野地空间分布

在荒野度分布数据的基础上, 将中国陆域空间划分为荒野地(荒野度大于等于3)、荒野缓冲区(荒野度等于2)、非荒野缓冲区(荒野度等于1)和人类长期活动区四种类型(图 6)。荒野地占全国陆地面积的22.47%, 主要分布在青藏高原、柴达木盆地、准噶尔盆地、塔里木盆地、内蒙古高原西部、大兴安岭和小兴安岭等区域。荒野缓冲区属荒野度较高但仍受一定程度的人类活动影响, 且经过生态保护最可能重新转变为荒野地的区域。非荒野缓冲区靠近人类长期活动区, 受人类影响较大, 如不加以限制和保护, 很可能转变为人类长期活动区。

自然保护地是在政策和法律支持下划定的自然生态和生物多样性保护区, 而荒野地是自然界留存的生物“最后的避难所”^[2], 自然保护地和荒野地空间分布关系的研究, 或可为自然保护地的自然环境现状分析, 和中国荒野地保护的必要性研究提供参考。本文参照戴文昱等^[51]的工作, 选取主要的自然保护地分布数据, 包括国家公园和自然保护区。由图 6, 位于藏北高原和罗布泊的自然保护区与荒野地有很高的重合度, 大小兴安岭、横断山脉和腾格里沙漠等地区的荒野地区域内, 涵盖众多面积较小的自然保护地。另外, 数量众多的小面积保护地分布在中国东南部的非荒野缓冲区, 此类地区距离人口聚居区比较近, 受人类活动影响较大, 因此是目前需要加强保护的自然保护地。

2.3 荒野地识别结果对比

以Mu等^[27]基于人类足迹法识别的荒野地为对比数据, 首先将此数据由1 km分辨率重采样为30 m, 另按照其划分类型, 将本文结果重新划分为荒野地(荒野度大于等于3)、未受破坏区(荒野度为2)和人类高度开发区(人类长期活动区和荒野度等于1的区域)。利用Combine函数运算, 提取两种数据三个类型的组合数, 并将九种组合分为四类, 结果如图 7。两种方法识别结果完全相同的区域面积占比超过80%, 仅本文方法或者对比数据单方面识别为荒野地的区域仅9%, 其中仅本文识别的荒野地主要分布在内蒙古高原、藏北高原南部、准噶尔盆地和横断山脉等地区, 仅对比数据识别的荒野地主要分布在柴达木盆地、巴颜克拉山地区、杭州湾、太湖及长江入海口区域。指标体系及评价方法的差异是导致识别结果不同的重要原因, 例如, 大型牧场多位于偏远地区, 荒野度比较高, 但本文引入牧场分布数据, 利用布尔运算判定其并非荒野地, 如柴达木盆地和巴颜喀拉山部分地区; 建立可达性指标可表征人类活动非线性扩散的规律, 作为乘积因子直接影响荒野度范围, 距离人类居住区近且易达的地区荒野程度更低, 如杭州湾和太湖等区域, 而距离居住区较远且极难到达的地区荒野度比较高, 如藏北高原南部地区等。

3 讨论与结论

构建科学合理的指标体系是准确识别荒野地空间范围的关键。以往的荒野连续谱法和人类足迹法建立一系列表征自然程度或人类压力的指标, 通过指标叠加量化得到可表征荒野程度的连续变量, 选取更多的指标类型有利于完善模拟结果, 其阈值的选择可直接影响荒野地空间范围的确定。马力等^[2]建立的决定—影响耦合关系模型, 从人类是否宜居出发区分荒野地和非荒野地, 该模型能有效规避人为确定阈值可能带来的风险。Ma等^[13]采用了可直接反应人类活动密度和强度的腾讯LBS(Location Based Services)位置大数据模拟荒野地分布, 此方法的精度更多依赖于数据源的质量, 但可避免建立模型时可能引入的误差。

本文从人类活动对自然影响的角度出发构建指标体系, 基于人类非持续性活动强度评估荒野度, 综合人类长期活动区界定可能的荒野区, 实现对中国陆域荒野地的识别。引入光污染等新型数据源, 选取较高精度的地理空间数据, 探索提高大尺度荒野地的识别精度。另外, 聚焦人类与自然的交互行为并建立了人类非持续性活动计算模型, 以评估荒野地质量。基于本文方法的评估结果与基于人类足迹法的荒野地数据进行对比, 在总体上具有较高的重合度, 并且在部分区域内表现出较优的识别能力。然而, 本研究也存在以下不足: 本文虽然选择了精度较高的DEM、土地覆盖和居民区等数据, 但牲畜分布和光污染等数据的分辨率并不高, 尽管采用了降尺度方法进行了数据处理, 但此结果并不能完全替代实际的30 m分辨率数据, 因此需要继续寻找更加适用的高精度数据集; 牧场是本模型用于识别非荒野地的一项重要指标, 本文采用包含牛、山羊和绵羊三种牲畜的牲畜分布数据, 并以草原为主要地表景观类型, 识别了草原牧场, 但由于牧场数据生产的局限性, 目前没有适用性较好的数据源, 因此忽略了例如马、骆驼等牲畜类型, 以及戈壁等牧场类型, 下一步工作中将关注更多数据源, 优化牧场识别模型以提高识别质量。