新时代背景下, 生态安全已经成为关系人民福祉和民族未来的大事, 生态安全格局的构建对于缓解生态保护和经济发展的矛盾具有重大的意义^[1-2]。生态安全具体指人类的生产生活和对土地资源的需求以及生态环境改变的适应方式等不会产生威胁的状态^[3-5]。土地利用生态安全是生态安全内涵的延伸, 指通过优化土地利用结构从而降低对生态环境变化和安全产生的影响, 主要方向以通过改变土地自然要素, 如土壤、水文等来改善土地生态系统的安全状态^[6-7]。生态安全格局是在人为技术手段和规划设计的前提下, 实现多层次优化方案, 以求实现最大的综合效益, 改善土地功能, 达到可持续发展的空间优化格局^[8-9]。Seppelt等^[10]在GIS平台和景观空间分异模型的前提下以Hunting Creek流域为研究区, 通过规划土地空间优化方案以及规范化肥的使用标准来改善土地的污染问题; Allan等^[11]在GIS平台下为保障小流域的水环境, 对小流域进行了空间优化。俞孔坚^[12]提出了适合国内的“景观安全格局”理论, 主要指生物资源迁徙所克服的最小景观阻力所做的功, 生态安全格局是将点、线、面相互叠加所组成; 陈利顶、傅伯杰等^[13-14]从“源”和“汇”角度将景观格局与土地生态过程结合在一起进行研究; 魏伟等^[15]采用最小累计阻力模型和预测模型构建了石羊河流域土地利用生态安全格局, 将生态安全格局与土地利用预测情景结合, 以期得到最优的土地利用生态安全格局。

《河北雄安新区规划纲要》中提到要坚持生态保护优先, 打造“太行山脉—渤海湾”与“京南生态绿楔－拒马河－白洋淀”生态廊道, 形成连山通海、南北交融的区域生态安全格局。鉴于此, 本文选取“太行山脉—渤海湾”与“京南生态绿楔－拒马河－白洋淀”两生态廊道上15个县(市)作为河北省典型样带, 利用最小阻力模型(MCR), 构建土地利用生态安全格局, 以期能够缓解该研究区生态保护和经济发展之间的矛盾, 为“京津冀”协调发展和“雄安新区千年计划”建设提供借鉴。

1 河北典型样带选取原则及依据 1.1 选取原则

参照IGBP/GCTE对陆地样带^[16-17]设置的原则, 本文在河北省布设典型样带, 具体位置为张家口蔚县, 保定市涞源县、易县、定兴县、徐水县, 雄安新区的容城、雄县、安新, 廊坊市的文安县和大城县, 沧州市的任丘市、青县、沧县、黄骅市和沧州市城区(图 1)。在选取样带时, 遵循以下原则: 1)典型性:在样带选取地域土地类型、地形地貌、土壤类型、气候(温度、降水等)等土地自然要素须具有典型梯度地带性和空间多样性, 即样带能够在最短的距离内穿过最全的生态景观类型且具有明显梯度地带性, 能够代表一定的区域。2)易用性:样带的选择过程符合实际调查结果, 具有现实意义, 即后期研究中的数据获取以及野外调查时容易实现, 并且对于热点区域、敏感地带、生态安全等的研究具有深远影响。3)完整性:样带的选取尽可能涵盖整个区域。比如:中国东北样带的选取是以研究全球气候变暖为目的选取的典型样带, 具有典型的气候、植被、土壤和土地利用过渡特征^[18-20]; 中国东部南北样带以热量与水热组合规律最为显著, 是全球独有的以热量为主要驱动力的全覆盖森林连续样带^[21-22]。

1.2 选取依据

河北典型样带位于114°13′~117°49′E, 38°05′~40°10′N, 起于西北部的蔚县山区, 止于东部的黄骅滨海地带, 直线距离300 km。地处河北腹地, 京津外围, 内含雄安新区, 既是粮食主产区、环京津重点水源涵养、生物多样性等重要功能保护区, 也是河湖滨海保护和盐渍化敏感地带。该样带具有典型的地貌、土地利用、气候、土壤等地带性特征。其中地形地貌由西北到东南主要跨越了间山盆地—山地—丘陵—山前倾斜平原—中部平原—滨海平原; 土地利用景观类型分布主要由地形因子决定了空间分布格局, 由西北向东南主要以草原、草甸、耕地—林地、灌木—农作物、乔木—盐荒草地、林果、农作物—盐荒草地、林果、农作物、湿地为主; 从气候要素看, 该样带属温带季风气候, 该样带由于经纬度跨越不明显, 所以气候差异变化在宏观尺度上相对不显著, 主要是以山地丘陵区与平原区过渡带为分界线, 山区年均气温小于4 ℃, 年累积降水量在430 mm以下, 平原区年均温高于10 ℃, 年累积降雨量高于500 mm; 土壤由西北向东南分为地带性与非地带性土壤, 其中地带性土壤分布为棕壤性土→棕壤→山地草甸土→栗褐土→淋溶褐土→褐土性土→石质土→粗骨土→石灰性褐土→潮褐土, 适宜林草地生长, 非地带性土壤为潮土、水稻土、灌淤土、沼泽土、盐化潮土等人为熟化的半水成土, 适宜农耕。另外, 该样带水系众多, 属海河流域大清河水系, 其主要为白洋淀及其支流, 以及蔚县的壶流河和中西部众多的人工水库、景观水面等。

2 数据来源和研究方法 2.1 遥感影像数据获取与处理 2.1.1 土地利用/覆被数据

本文采用的遥感数据均下载于地理空间数据云(http://www.gscloud.cn/), 包括5景图像(获取月份均为当年7—9月植被最为旺盛时段), 为保证遥感解译时空间尺度一致性和土地利用数据精度, 本文所用影像空间分辨率均为30 m。利用ENVI 5.3软件对影像进行波段组合、几何校正(误差值不超过0.5个像元)、图像镶嵌及行政区裁剪等前期基本处理, 并进行投影变换。根据土地覆被分类系统分为6类, 即耕地、草地、林地、湿地、人工表面、其他用地(图 2), 解译精度达90%以上。

2.1.2 NDVI数据

首先对Landsat影像进行数据预处理, 经过辐射定标、大气校正消除辐射误差, 后进行几何校正, 经过镶嵌裁剪得到预处理后的研究区遥感影像。

NDVI是表示植被生长好坏、覆盖度多少的指标, 值域为[-1, 1], 小于0代表云、水、雪等高反射的附着物, 0代表岩石或裸土, 大于0代表有植被覆盖。NDVI计算公式为:

$ \mathrm{NDVI}=\left(\rho_{\mathrm{NIR}}-\rho_{\mathrm{R}}\right) /\left(\rho_{\mathrm{NIR}}+\rho_{\mathrm{R}}\right) $

(1)

式中: ρ_NIR是近红外反射率; ρ_R是红光反射率, 在OLI影像中分别为波段5和波段4。

2.1.3 DEM数据

本文DEM数据下载于地理空间数据云(http://www.gscloud.cn/)的GDEMV2 30 m分辨率数字高程数据, 共有8幅影像, 轨道号分别为114/40、115/40、114/39、115/39、115/38、116/39、116/38、117/38。运用ENVI 5.3软件进行影像镶嵌及行政区裁剪等前期基本处理。

2.2 其他获取与处理

1) 分别收集到了研究区15个地级市(县)1:25万土壤扫描图, 运用MAPGIS、ARCGIS等软件经过数字化得到样带土壤图。2)高程、坡度等数据均通过DEM提取得到。3)水体图通过ARCGIS软件对土地利用覆被进行提取得到。4)收集到了各县级行政区的土地利用现状图和行政区界线(来源于全国第2次土地调查数据库), 以用于遥感数据解译和数据裁剪。5)气温数据来源于中国科学院资源环境科学数据中心(http://www.resdc.cn)。6)降雨数据来源于中国气象科学数据共享服务网(http://data.cma.cn/), 经过克里金插值得到。7)社会经济数据来源于《河北农村统计年鉴》及相关地市及统计年鉴。8)《河北雄安新区总体规划(2018—2035年)》《京津冀协同发展规划纲要》。以上数据在运用过程中均根据需要进行了相同的格式转换和投影转换, 栅格数据像元大小均重采样为30 m×30 m。

2.3 生态安全格局构建方法及步骤 2.3.1 “源地”提取

1) 生态用地保护“源”的提取

本文主要侧重于该样带土地利用模式对京津冀和雄安新区的发展起到生态保护作用, 尤其是连片的林地、草地和湿地等。该研究区成片林地主要集中分布在上游的蔚县、涞源和易县的山地丘陵区, 是国家未来重点关注的水源涵养林保护区, 因此作为生态保护的重点区域; 草地主要分布在蔚县山间盆地外围边缘地带、涞源东部中、低山地带以及易县山前倾斜平原地带, 也是作为生态恢复和防风固沙的主要区域; 此外, 白洋淀及其上游萍河、瀑河、漕河、唐河、府河、白沟引河等支流以及人工建设的水库、城市景观水面等对干旱内陆地区来讲是维系生态系统得以发展的“生命线”, 因此也是划定生态保护区时考虑的重点(图 3)。“生态源”是构成高安全水平的生态用地, 是保障自然生态系统的最小土地底线, 原则上任何城市扩张行为不得侵占这类生态用地。

2) 城镇用地扩张“源”的提取

城镇用地扩张的“源”是城镇扩张的核心动力。在生活用地中, 主要指各类人工表面用地, 本文通过河北典型样带2017年土地利用图提取各类人工表面用地得到河北典型样带城镇源地面积为1 650.04 km², 占研究区总面积的8.71%, 呈点状分布, 且平原地带居民点密集度高于山区(图 4)。

2.3.2 最小累积阻力面构建

从土地系统的生态过程角度来看, 生态和生活用地是一对具有竞争性质的关系, 突破这种竞争关系需要克服阻力, 而克服这种阻力在不同景观空间格局下所做的功不同, 所以二者在竞争过程中克服的阻力大小也不同, 为使研究具有意义, 必须建立相同标准下的阻力评价体系, 但所赋分值必须是相反的。

本文从地形地貌、土地景观、土壤、生态敏感性和生态系统服务价值5个方面建立阻力表面评价体系(表 1)。以系统论^[23]为基础, 将上述因子分为土地内在固有属性和外在生态属性, 其中土地固有属性指土地自然的本底要素与要素之间的相互作用, 决定了土地结构和功能的外部特征, 包括地形地貌、土壤, 反映土地利用空间格局分布受固有属性的约束; 土地外在生态属性是指土地系统功能和外部环境发生作用时所表现出来的特征, 包括土地景观类型、生态敏感性和生态系统服务价值, 反映该地区受生态约束的能力^[23-25]。

在构建过程源的阻力表面时, 把评价指标分为1、2、3、4、5等级, 其值越小, 代表所克服阻力越小, 反之越大, 赋分时除参考Li等^[26]对植被、土壤、景观类型的赋值方法外, 主要通过AHP法确定指标相对重要度。

1) 生态源阻力因子的选择及权重的确定

生态源阻力因子是从土地内在固有属性和外在生态属性两方面选取的, 其中内在属性包括高程和土壤类型, 地形用海拔分区来表达; 土壤类型依据土壤质地和肥力情况进行等级划分。外在生态属性包括土地景观类型、生态敏感性和生态系统服务价值。土地景观类型指上文解译的六大地类; 生态敏感性分级参考《HJ19—2011:环境影响评价技术导则:生态影响》所提供的分值和权重, 采用坡度、植被覆盖度、气温、降水4个要素进行评价分级; 在对生态系统服务价值进行评价时, 采用以NDVI和空间各像元距水体像元的距离为基准的分类矩阵法来分析^[27]。生态源阻力表面的评价指标体系如表 2所示。本文通过借鉴前人^[25]研究成果, 结合专家打分法和AHP法确定各指标权重, 将各个阻力表面空间分布图经过Arcgis软件的空间分析功能的栅格计算器进行求和后根据值域大小分为5个等级, 综合后的结果如图 5所示。

2) 城镇扩张源阻力因子的选择及权重的确定

对生活用地构建阻力表面评价指标时赋予分值与生态源地相反, 同时指标选取也通过内在固有属性和外在生态属性两方面选取。详细指标见表 3, 综合后的阻力值空间分布如图 6所示。

2.3.3 最小累积阻力模型

本文参考刘孝富等^[25]研究结果, 将研究区土地利用的优化划分为适宜生态的用地和适宜生活的用地两大类, 为平衡生态用地保护和经济发展之间的矛盾, 在相同标准下选取生态保护用地扩张和城镇用地扩张阻力因子, 并赋予阻力分值及权重, 通过MCR模型分别计算生态用地扩张和城镇用地扩张阻力面。

MCR模型由Knaapen等^[28]于1992年提出, 最初用于景观规划, 后由俞孔坚等^[29]用于生物保护的生态安全格局构建中。其基本原理是指生物资源从“源点”到“源地”克服景观阻力所做的功, 反映了生物资源迁徙过程的一种可达性^[22]。公式如下:

$ {\rm{MCR}} = f\min \sum\limits_{j = n}^{i = m} {{D_{ij}}} \times {R_i} $

(2)

式中: MCR为最小累积阻力值; f为最小累积阻力与生态过程的正相关函数; D_ij为从源点j到源地单元i的空间距离; R_i为单元i对某物种运动的阻力系数; $ \mathop \sum \limits_{j = n}^{i = m} $为源j与单元i之间所穿越所有单元的距离和阻力的累积; min为单元i对于不同的源取累积阻力最小值。

2.3.4 生态安全格局分区

本文参考刘孝富等^[25]的研究结果, 生态安全分区通过采用最小累积阻力差值与栅格面积关系曲线确定, 把频率发生急剧变化的点称为拐点, 反映其两侧的阻力值出现较大突然变化, 显示出高度的景观异质性, 将拐点处的阻力值视为临界值, 根据阻力值频率变化特点与实地情况分析划分生态安全区。本文运用Arcgis软件平台对“生态扩张阻力面”和“城镇扩张阻力面”进行差值计算, 以二者差值进行生态安全分区。公式如下:

$ {\rm MCR}_{差值}={\rm MCR}_{生态用地扩张}-{\rm MCR}_{城镇用地扩张} $

(3)

式中: MCR_{生态用地扩张}为生态扩张综合最小累积阻力值, MCR_{城镇用地扩张}为城镇用地扩张综合最小累积阻力值。当MCR_差值 < 0时, 即该区域生态用地扩张阻力小于城镇用地扩张阻力, 适合生态用地扩张, 并且差值越小越适宜生态用地扩张; 当MCR_差值 > 0时, 该区域生态用地扩张阻力大于城镇用地扩张阻力, 适合城镇用地扩张, 并且值越大越适宜城镇用地扩张; 当MCR_差值=0时, 为该区域的生态用地扩张与城镇用地扩张的分界线。

3 结果与分析 3.1 最小累计阻力面的生成

根据公式(1), 分别计算生态扩张最小累积阻力面和城镇扩张最小累积阻力面。由图 7可知, 生态源与城镇扩张源的最小累积阻力空间分布差异较为明显。生态源最小和最高累积阻力值分别为0和29 708.5。其中高阻力值区主要分布在各县域城镇聚集区及乡镇、农村等人为活动强烈地区, 且这类区域基本处于海拔低、气温高、降雨量大, 经济发展和城镇扩张较快, 以农业生产和城镇建设为主的平原地区, 生态保护过程的阻力较强, 其内部生态流和能量流传递十分有限, 因此其累积阻力值较高, 是生态治理的重点地区。低阻力值在整个样带均有分布, 主要分布在蔚县、涞源、易县的西北部山地丘陵区, 以及中部平原白洋淀流域外围和各个湿地外围缓冲地带。生态源低阻力值分布区域基本涵盖两个特点:一是河网密布, 水源较为丰富, 尤其在山前倾斜平原和中部平原地带, 上游出山口地区水量较大并汇聚于白洋淀, 地处国家自然保护区的核心区; 二是这些地方植被覆盖度较高, 上游山地丘陵区是典型的水源涵养林分布区, 中下游平原地带是重要的农耕区, 白洋淀流域及各大支流从西北至东汇聚, 沿河林、草、耕地分布较多。以上两个特点使得河北典型样带内部生态流和能量流交换传递的速度加快, 因此其生态源最小累积阻力较低, 是生态保护的重点区域。

城镇扩张源的累积阻力最低值和最高值分别为0和31 831.0。低值区主要分布在蔚县山间盆地和平原地区, 这些地区既有丰富的水源, 又有便利的交通, 也有较为发达的农业系统和经济繁荣的城镇, 承载着河北典型样带绝大多数人口, 因此对于城镇发展而言, 所受阻力值较小, 适合城乡建设用地优化布局; 而高值区主要分布在海拔相对较高的蔚县、涞源、易县的山地丘陵区, 以及中部平原白洋淀流域和黄花滨海地带, 其原因是山地丘陵区海拔高、气温低、植被茂密, 为人居环境不适宜区, 其他累积阻力高值区生态环境恶劣, 交通水源不便, 也不适合人类居住, 因此这些地方对于城镇发展和城镇扩张源而言, 其内部生态流、能力流所受阻力较大, 其最小累积阻力也较高。

从生态源和城镇扩张源值域分布来看, 生态源整体累积阻力小于生活扩张源累积阻力, 说明在河北典型样带城乡建设用地生态流和能量流内部传递交换的阻力远大于生态用地。综合累积阻力高值区分布范围来看, 无论是生态源还是扩张源, 在雄安新区以及上游定兴和徐水地区, 其值都很高, 证明这些地方生态阻力很强, 应是生态治理的重点地区, 从流域可持续发展来说, 生态保护重要性大于城镇建设, 该样带应加强生态保护。

3.2 生态安全格局构建 3.2.1 生态安全分区阈值确定

依据公式(2), 运用Arc GIS中的栅格计算器将生态扩张和城镇扩张的累积阻力表面进行差值计算, 最终生成最小累积阻力差值表面, 如图 8所示。从图中可以看出, 生成的最小累积阻力差值范围为-31 138.74~27 695.06, 累积阻力差值大于0的区域主要分布在建设用地集聚、人口密度较大、水源较丰富的区域, 其中平原属于高值密集区, 该区域人类活动剧烈, 对土地的利用程度较强, 对于生态扩张阻力较大; 累积阻力差值小于0的区域, 主要分布在蔚县东南部、涞源山地区域以及易县东部山地丘陵区, 该区域海拔较高、气温低、坡度大, 土壤主要为山地草甸土、棕壤、栗褐土、褐土性土等不宜耕种的土壤, 不适于人类居住和建设用地的开发, 因此是主要生态源地所处的区域。

根据图 8, 以最小累积阻力值差值与栅格面积的变化的突变点来确定阈值并划分土地利用生态安全分区。在前文叙述基础上, 将生态适宜用地进一步划分为生态核心区和生态缓冲区, 适宜建设用地进一步划分为生态过渡区、生态可占区和生态治理区。在累积阻力差值小于0的部分, 于MCR_差值等于-10 225.163 2处斜率发生突变; 在累积阻力差值大于0的部分, 于MCR_差值等于1 725.453 4和6 781.483 5处斜率发生突变。突变点前后的土地变化较大, 因此, 将MCR_差值等于-10 225.163 2处、MCR差值等于1 725.453 4和6 781.483 5处作为分区点, 划分结果见表 4。

3.2.2 土地利用生态安全分区

根据分区阈值划分标准得到生态安全格局分区的数值及空间分布, 如表 5和图 9所示。其中, 生态核心区主要分布在样带西北部, 具体位于蔚县东南部2 000 m以上、涞源和易县1 500 m以上的高海拔区域, 面积约1 082.42 km², 占比5.77%。这里森林密布、气温低、坡度大、河网沟壑密集, 是河北典型样带自然保护区核心地带, 严禁从事矿产开发和破坏环境等行为。

生态缓冲区主要分布在样带中部和东部, 具体位于蔚县北部、东南部, 涞源全域和易县中东部的中低山地区, 以及国家自然湿地保护区白洋淀和黄骅滨海地带, 总面积最大为6 943.63 km², 占比37.01%, 这里主要由灌木林、天然草地以及湿地组成, 是连接上、中、下游的媒介, 也是保障整个河北样带生态安全的屏障, 因此应该在保护的基础上合理利用, 防止出现过度开发。

生态过渡区主要分布在样带中部和东南部, 具体位于蔚县间山盆地、涞源县城和易县东部山前倾斜平原区, 以及中部平原和滨海平原大部分地区, 总面积达3 189.25 km², 占比17%。该区由部分生态型用地和一般型生态用地以及大部分耕地组成, 这里生态保护、农业生产和城镇建设矛盾较为突出, 也是生态扩张与城镇扩张的交界处。

生态可占区主要分布在样带中部和东南部, 具体位于生态过渡区和生态治理区的过度地带, 主要包括现已开发的城市建成区和重点建设的乡镇以及周边具有开发和发展潜力的区域, 是今后作为城镇扩张重点开发建设区域, 应合理规划。同时该区外围分布着部分质量较差的耕地, 也包括大部分低效建设用地, 也是今后进行土地整治和开发的重点后备区域, 面积为5 835.34 km², 占比31.10%。

生态治理区主要分布在样带中部和东南部, 具体位于白洋淀上游的徐水县东北部、定兴县东南部、容城和雄县地区, 沧县东部以及零散分布在平原的各个生态可占区的外围, 这以耕地和湿地为主, 不但是生态脆弱区, 也是未来进行生态恢复、生态产业发展的重点地区, 面积为1 710.70 km², 占比9.12%。该区应加强治理力度, 积极探索自然生态恢复, 并尝试生态产业和产业链发展。

4 讨论与结论

彭建等^[4]的研究认为, 生态安全格局作为沟通生态系统服务和人类社会发展的桥梁, 目前被视为区域生态安全保障和人类福祉提升的关键环节。本文以描述生态过程的最小累积阻力模型构建了以最小累积阻力差值表面进行生态安全分区模式, 利用此种方法研究生态安全格局也已较为成熟。刘孝富等^[25]、李颖丽^[3]的研究结果, 以最小累积阻力值差值与栅格面积的变化的突变点来确定阈值并划分土地利用生态安全分区, 将河北省典型样带划分为5个区域, 分别为生态核心区、生态缓冲区、生态过渡区、生态可占区、生态治理区。因此, 在构建生态安全格局的基础上, 针对样带不同的生态安全区问题提出可行性的生态环境安全可持续发展对策。生态核心区和生态缓冲区共同构成中水平的生态安全格局, 因地制宜大规模植树造林, 尤其是在裸岩、裸土等区域, 提高植被覆盖度, 减少自然灾害, 严格禁止任何破坏行为; 生态过渡区以生态保护和农业生产为主, 属于限制建设区, 与生态核心区和生态缓冲区共同构成高水平生态安全格局, 在保障水土资源基础上进行农业经济生产开发, 严格控制向生态用地扩张。生态可占区是该条带的核心建设与经济发展区域, 以生态建设与经济发展并重, 协同雄安新区的发展规划, 根据城市人口、经济等需求进行因地制宜开发和合力建设, 打造绿色智慧城市; 生态治理区是该条带亟待解决的重点区域, 恢复淀泊水面, 实施退耕还淀, 增强白洋淀生态自我修复能力, 推动区域环境协同治理, 严守土壤环境安全底线, 优化生态安全屏障体系, 提升生态系统质量。

本文以土地利用现状确定生态源地和城镇扩张园地, 避免主观人为和模型误差的影响; 在指标选取上以系统论为基础, 从土地内在固有属性和外在生态属性两方面进行选取, 结合MCR模型, 进一步加强了“源”与生态过程的联系。但本文也存在着不足, “源”的位置、空间距离、阻力面等的影响因素很多(人文驱动、社会经济等); 另外在赋指标分值和指标权重测算过程中, 人为主观影响较大, 对生成的最小累积阻力表面有一定的误差; 数据来源全部为遥感数据, 没有外业调查和对指标进行验证, 对结果有一定的影响。综上部分不足应有待进一步研究。