20世纪90年代以来, 随着城市化进程不断加快, 资源环境与土地利用的矛盾逐渐尖锐, 生态与环境问题日益凸显, 科学开展生态风险评估并识别不同区域风险等级, 实现精准施策, 成为了亟需解决的问题^[1—2]。生态风险评价研究最早始于美国20世纪70年代, 1992年美国环境保护署将其定义为由一种或者多种因子导致的可能或者正在发生的不利生态影响的过程^[3]。国内外学者针对生态风险评价开展了大量的研究工作, 生态风险评价可从多角度以及多尺度开展, 包括水体中的微生物生态风险^[4]、沉积物中的重金属生态风险^[5]、碳排放风险^[6]、土地整理后的生态风险^[7]以及城市和流域生态风险^[8—9]等。在流域尺度的生态风险评价手段上, 主要的评价方法有最小阻力模型^[10]、景观生态风险评价模型^[11]、基于压力-状态-响应模型^[12]、综合生态敏感性指数^[13]等, 用于指导流域空间管控及分区划定^[14]。景观生态风险评价是利用景观格局指数, 建立景观结构与生态过程相关关系, 反映景观格局中的生态风险。从景观生态风险研究对象上来看, 目前以流域单元为研究对象的较多, 国内外学者对研究区域开展景观生态风险评价研究, 划分生态风险网格, 通过数学模型计算得出景观生态风险的高中低分布, 但进一步结合研究成果特别是识别高风险区域后再进行局部研究的较少。密云水库是首都最重要的地表饮用水水源地和水资源战略储备基地, 对保障首都供水安全具有重要意义^[15], 同时密云水库地处北京市生态涵养区, 流域内遍布自然保护区和生态保护红线, 是重要的生态屏障。习近平总书记在给建设和守护密云水库乡亲们的回信中强调, 密云水库作为北京重要的地表饮用水水源地、水资源战略储备基地, 已成为无价之宝^[16]。2021年中共北京市委生态文明建设委员会印发的《北京市密云水库流域水生态保护与发展规划(2021—2035年)》构建了北京市密云水库流域“一库一环三区多廊”的水生态空间保护格局, 近日北京市规自委编制的《北京市密云水库上游地区空间保护规划(2020—2035年)》, 构建了“1+2+2+179”整体空间结构, 统领水库上游地区形成与生态资源环境相适应的空间保护和发展格局。2016年, 北京市划定了密云水库饮用水水源保护区范围, 将密云水库流域划分为一级保护区、二级保护区以及准保护区, 其中一级保护区为密云水库环库公路以内(荞麦峪西侧至口门子村、城子以南至黄土洼以北、前保峪岭至老爷庙背水一侧及鲶鱼沟南背水一侧划定的区域除外), 以及环库公路以外的进水地带, 二级保护区为一级保护区之外至密云水库的向水坡范围以内以及密云水库调节池的汇水范围以内, 准保护区为二级保护区以外上游河道的流域。从多时空、多层级角度开展密云水库流域生态风险评价, 识别流域不同时期生态风险演变情况, 对流域可持续发展和未来土地规划具有重要意义。

本研究从景观角度, 基于景观干扰度、景观脆弱和景观损失度构建1990—2018年4个时期生态风险评价模型, 从流域尺度、水源保护区尺度, 划分流域生态风险等级, 同时基于莫兰指数理论, 分析空间自相关关系, 支撑以网格为单元的精细化管理, 为流域生态环境保护和修复更具生态价值和战略意义提供借鉴。

1 材料与方法 1.1 研究区域概况

密云水库流域面积1.58万km²(图 1), 其中占河北省面积1.23万km², 占北京市面积0.35万km², 河北省涉及张家口市沽源、赤城, 承德市丰宁、兴隆、滦平5个县、北京市涉及延庆、怀柔、密云3个区, 总人口183.5万人。密云水库流域属北京市五大水系中潮白河水系, 有潮河、白河、汤河、天河、安达木河等主要河流。流域属暖温带季风型大陆性半湿润半干旱气候, 全年平均气温为10.8℃, 主要土壤类型为褐土和棕壤土, 降水分布一般是从东南向西北递减, 主要降雨集中在6—8月汛期, 降水量300—700mm^[17]。

1.2 数据来源

为探究密云水库流域生态风险时空变化趋势, 选择研究时段为1990—2018年。数据来源为中国科学院资源环境数据中心(https://www.resdc.cn), 包括1990、2000、2010和2018年全国土地利用遥感监测数据(栅格类型, 分辨率30m), 该数据是基于美国陆地卫星Landsat TM影像, 通过人工目视解译生成, 土地利用类型包括6个一级类型以及25个二级类型。根据研究需要, 采用ArcGIS 10.2软件对一级类进行分析, 即耕地、林地、草地、水域、建设用地和未利用地六种用地类型(图 2)。

1.3 研究方法 1.3.1 景观格局指数计算

景观指数能够高度浓缩景观格局信息, 并反映其结构组成和空间配置某些特征, 是景观异质性的具体表现, 通过景观格局指数计算, 可以反映景观大小、形状、排列及生态过程等, 一般分为三种尺度, 即斑块、类型和景观尺度。根据前人对密云水库流域以及水质和景观格局的相关研究^[18—20], 用Fragstats 4.2软件进行景观格局指数计算, 选择对水质具有显著影响的景观格局指数, 在景观尺度下选择形状指数(LSI)、蔓延度指数(CONTAG)、分离度指数(SPLIT)和香农多样性指数(SHDI), 在类型尺度下选择斑块数(NP)、斑块密度(PD)、最大斑块指数(LPI)(表 1)。

1.3.2 生态风险评价模型构建

基于景观生态学理论, 选取能够反映区域空间生态风险的相关指数, 通过计算生态风险单元景观干扰度、脆弱度、损失度等生态指数, 构建生态风险评价模型, 揭示流域内生态风险变化趋势。

(1) 生态风险单元划定

将流域划分等间距生态风险网格, 细化生态风险分布规律和差异, 流域内平均斑块面积约为1.81km², 按照平均斑块面积的5倍创建渔网^[21], 即3km×3km网格, 将整个流域划分为1867个风险单元。

(2) 景观干扰度指数(S_i)

S_i反映不同景观类型受外部干扰程度, 通过计算景观破碎度指数、景观分离度指数和景观优势度指数赋权计算得出。表达式如下：

式中Ci为景观破碎度指数, n_i为景观类型i的斑块数, A_i为景观类型i的面积, A为所有景观的面积, Ni为景观分离度指数, Di为景观优势度, Q_i为景观类型i出现的样方数与总样方数的比值, M_i为景观类型i的斑块数与总斑块数的比值, L_i为景观类型i的斑块面积与样方总面积的比值, a、b、c为权重值, a+b+c=1根据相关研究, 分别赋值0.5、0.3、0.2。

3) 景观脆弱度指数(F_i)

F_i指景观受外部干扰的脆弱程度, 根据相关研究^[22], 对各类景观脆弱程度进行赋值为：耕地4, 林地2, 草地3, 水域5, 建设用地1, 未利用地6, 赋值后进行归一化处理。

4) 景观损失度指数(R_i)

R_i可通过每一种景观的干扰度和脆弱度相乘计算得出。表达式如下：

S_i为景观干扰度指数, F_i为景观脆弱度指数。

5) 生态风险评价模型(ERI)

根据景观干扰度、景观脆弱度、景观损失度计算生态风险指数, 将景观空间结构转化为空间生态风险, 即构建生态风险评价模型。表达式如下：

ERI_i为生态风险指数, n为生态风险单元内景观类型数量, A_i为第i类生态风险单元景观面积, A为景观总面积。

1.3.3 空间相关性分析

通过计算Moran′s I分析流域全局空间相关性, 从而揭示流域内聚集和离散程度^[23]。表达式如下：

式中Y_i和Y_j为变量在相邻配对空间单元的变量, ω_ij为空间权重矩阵, Y为平均值。Moran′s I＞0表示空间正相关性, 其值越大, 空间相关性越明显, Moran′s I＜0表示空间负相关性, 其值越小, 空间差异越大, 否则, Moran′s I = 0, 空间呈随机性。

局部Moran′s I指数可以反映区域与相邻区域之间的相关性及显著程度, 局部空间自相关指标(LISA)聚类地图可以可视化分析相邻区域间差异程度^[24]。表达式如下：

式中, n为样本数量即生态风险单元的个数, S²为空间变量的统计方差。当I_i＞0时, 表示局部相邻空间正相关聚集(“高-高”、“低-低”), 当I_i＜0时, 表示局部相邻空间负相关聚集(“高-低”、“低-高”), 当I_i=0时, 表示相邻空间不显著, 呈随机分布。

2 结果与分析 2.1 景观时空变化 2.1.1 土地利用变化

根据四期遥感图解译的土地利用数据, 分析1990—2018年四期土地利用面积和变化情况。如表 2所示, 密云水库流域内主要土地利用类型为林地、草地和耕地, 三者约占总面积的97%, 其中林地面积占比最大。从时间变化来看, 耕地面积呈现逐年减少的趋势, 1990—2018年耕地减少了3.77%, 这与国家号召开展的“退耕还林”政策及北京市补偿河北省开展的“稻改旱”工程相关^[25], 一定程度上调整了土地结构和农业种植结构；水域用地呈现先增加后减少的趋势, 在2000年时水域面积最大为310.93km², 2010年水域面积比2000年减少了85.03km², 这与密云水库流域年降雨量及潮河、白河径流量和入库当年流量有关, 密云水库入库径流量在2000后表现出显著性减小^[26], 2018年比2010年增加了29.58km², 这与2014年后南水北调江水进京关系密切；随着城市建设的扩张, 未利用地逐年减少, 建设用地的面积逐年增加, 2018年建设用地的面积是1990年建设用地面积的3.5倍, 这是城镇经济社会发展的必然趋势。

2.1.2 景观指数变化

如表 3所示, 密云水库流域在1990—2018年间, 斑块数(NP)和斑块密度(PD)呈逐渐上升趋势, 分别增加了6.46%和6.78%, 可以推断随着流域内受人为干扰和城市建设扩张, 景观破碎化日益严重, 各类景观较为分散^[27]；最大斑块指数(LPI)呈现先增加后降低的趋势, 说明优势景观斑块控制作用在1990—2010年逐步增强, 在2010—2018年间随着景观破碎化严重, 优势斑块作用减弱；景观形状指数(LSI)呈现逐渐上升趋势, 增加了2.12%, 说明景观的形状变得越来越不规则和复杂；分离度指数(SPLIT)呈现先降低后增加的趋势, 特别是在2018年最大为47.78, 说明在2018年景观分离程度最强, 受人为干扰的破碎度最强烈；蔓延度指数(CONTAG)呈逐渐下降趋势, 减少了1.77%, 说明景观的团聚程度和延展趋势逐年降低, 也反映了破碎化程度加大；香农多样性指数(SHDI)呈逐渐上升趋势, 说明景观多样性增加。综上, 从景观格局指数分析, 流域受人为干扰影响逐年增加, 景观格局趋于复杂和分散, 破碎化加剧。

2.2 生态风险评价 2.2.1 流域生态风险时空变化

运用ArcGIS提取每个生态风险网格用地类型面积并用生态风险评价模型进行风险值计算, 运用地统计分析模块进行函数拟合和方差分析, 确定选用普通克里金球状模型进行拟合插值得到1990—2018四个时期生态风险时空分布, 综合考虑生态风险值可对比性, 基于自然断点法, 统一将生态风险分为5个等级对应的生态风险区, 即低风险区(ERI≤0.036)、较低风险区(0.036＜ERI≤0.043)、中风险区(0.043＜ERI≤0.056)、较高风险区(0.056＜ERI≤0.084)和高风险区(ERI>0.084), 得到生态风险时空分布图(图 3)并统计各风险等级面积及占比(表 4)。在时间尺度上, 1990—2018年生态风险值的最大值呈逐年上升趋势, 由1990年的0.23升高至2018年的0.36, 但低风险区和较低风险区的面积有所增加, 2018年低风险区和较低风险区面积总和比1990年增加了2255.35km², 占比14.28%, 且高风险区域面积减少了163.38km², 占比1.03%, 说明生态风险呈减弱趋势, 生态保护成效显著, 仅有局部地区生态风险值较高。在空间分布上, 1990年和2000年生态风险区域分布相似, 较高风险区主要分布在研究区北部(丰宁县、沽源县), 西南部(赤城县), 高风险区分布在研究区南部密云水库周边(密云区、兴隆县), 2010年和2018年生态风险分布相似, 主要以中低风险为主, 仅在研究区南部边缘出现小面积高风险区域。通过分析不同地类在各风险等级面积分布情况(表 5), 四个时期, 耕地、林地、建设用地和未利用地主要占比相似, 其中, 建设用地和未利用地主要分布在较低风险区和中风险区, 林地主要分布在较低风险区和低风险区；在1990年和2000年, 草地和水域用地主要分布在较低风险区和中风险区, 而在2010年和2018年, 二者主要分布在较低风险区和低风险区, 中风险占比减少, 生态风险逐渐减弱。整体上从用地占比来看, 高风险和较高风险占比相对较多的是耕地、水域和建设用地, 但也呈现逐渐减少趋势。

2.2.2 水源保护区生态风险时空变化

通过对流域生态风险时空变化分析得出中高风险主要分布在密云水库周边, 因此开展水源保护区生态风险评价是十分必要的。运用ArcGIS按掩膜提取功能, 分别提取一级、二级水源保护区生态风险栅格数据, 并运用自然断点法进行分级, 分别得到低风险区(一级保护区：ERI≤0.038, 二级保护区：ERI≤0.032)、较低风险区(一级保护区：0.038＜ERI≤0.041, 二级保护区：0.032＜ERI≤0.037)、中风险区(一级保护区：0.041＜ERI≤0.043, 二级保护区：0.037＜ERI≤0.042)、较高风险区(一级保护区：0.043＜ERI≤0.046, 二级保护区：0.042＜ERI≤0.045)和高风险区(一级保护区：ERI>0.046, 二级保护区：ERI>0.045)(图 4、图 5)。如表 6所示, 一级水源保护区高风险面积在1990年最大, 面积为149.65 km², 占一级水源保护区面积的54.82%, 且集中分布水源保护区中心位置, 在2018年最小, 面积为55.20 km², 占比20.22%, 空间分布上由中心向东北角转移, 低风险区面积逐渐增大, 在2018年最大, 面积为121.19 km², 占比44.39%。二级水源保护区高风险区在1990年最大, 面积为72.57 km², 占二级水源保护区面积的21.34%, 在2018年最小, 面积为32.70 km², 占比9.62%, 空间上与一级水源保护区高风险区集中连片；低风险面积逐渐增大, 面积由1990年的33.15 km²增加至2018年的86.46 km², 比例由9.75%增加至25.43%。整体上一二级水源保护区高风险区分布较为集中, 且高风险区面积逐渐减少, 生态风险趋于减弱, 与全流域生态风险变化趋势相同。

2.2.3 流域生态风险转移特征

运用ArcGIS对1990年和2018年土地利用数据进行空间叠加, 分析密云水库流域生态风险区转移情况。如表 7所示, 全流域生态风险区面积转移6346.8 km², 占比40.17%, 其中转出面积最大的是中风险区, 转出2620.53 km², 占流域面积的16.59%, 占转出面积的41.29%, 转入面积较大的是较低风险区, 转入3493.42 km², 占流域面积的22.11%, 占转入面积的55.04%。生态风险区之间转移最大的分别是中风险区转为较低风险区(2337.82km²)和较低风险区转为低风险区(1292.11 km²)；生态风险转移比例最大的分别是较高风险区向中风险区转移76.30%和中风险区向较低风险区转移51.21%, 进一步说明生态风险呈现下降趋势, 生态安全有所提升。

2.3 生态风险空间关联

运用GeoDa空间自相关软件构建以共边或共点为邻接的皇后权重矩阵(Queen), 分析1990—2018四个时期生态风险小区空间自相关情况, 计算Moran′s I指数(表 8)并绘制LISA聚类地图(图 6)。通过999次置换检验(P=0.001), Moran′s I大于0, 说明生态风险指数在空间上呈现正向的相似关系, 四期指数值先上升后下降, 在1990年, Moran′s I最大为0.322, 说明在1990年空间聚集程度最强, 随着社会经济发展和土地利用变化, 生态风险的聚集性减弱, 逐步趋于离散。通过Lisa聚类地图分析, 研究区域除不显著单元外, 正相关类型(“高-高”、“低-低”聚集单元)占比18.80%—22.87%, 负相关类型(“高-低”、“低-高”聚集单元)占比1.61%—2.14%, 整体上以“冷点”区域(“低-低”聚集单元)为主, 在2000年和2018年“冷点”区域个数最多, 分别为279和288个, 占比14.94%和15.43%, 且集中分布在北京市境内(延庆区北部、怀柔区南部、密云区西部)、零星分布在河北省境内(赤城县、丰宁县、兴隆县)；“热点”区域(“高-高”聚集单元)呈现逐年增加趋势, 四个时期“热点”区域主要出现在密云区南部、兴隆县南部、丰宁县北部, 且多数位于流域边缘, 在2018年最大, 为139个, 比2010年增加31.13%, 体现在赤城县西南部集中出现了“热点”区域。

3 讨论

在时间层面, 尽管随着人为干扰和城市建设扩张, 景观呈现了逐步破碎的趋势, 但景观多样性也逐年增加；研究区域位于水源保护区, 主要土地利用类型为林地和草地, 从“源汇”景观理论分析, 林草地对水源保护和生物多样性的生态过程演变起到了正向促进作用, 是“源”景观, 同时随着密云水库战略水源地的重要性逐年提升, 无论是“稻改旱”工程、“南水北调”工程、“百万亩造林”工程还是各项水源保护政策向该区域倾斜, 都有助于该流域生态逐渐向好, 高风险区域逐渐向低风险区转移, 生态安全逐渐提高^[28], 因此区别于其他流域^[29]由于城镇大面积扩张和粗放管理带来的生态风险加剧。在空间层面, 景观生态风险空间差异明显, 在流域尺度呈现了“边缘高、中间低”的分布规律, 而在水源保护区尺度上, 呈现了“中间高、边缘低”的分布规律, 其主要原因为流域尺度上流域边缘主要景观类型为水域和耕地, 虽然二者景观干扰度指数较低, 但景观脆弱度指数较高导致生态损失度相对较高, 导致生态风险值较高, 而水源保护区尺度中间景观类型为水域, 边缘为林地, 林地的生态风险相对较低, 且在不同时期水面面积变化, 高生态风险区也随之变化。生态风险分布和Lisa聚类地图可以看出, 流域变化趋势相近, 高风险区和高值都出现在相同区域, 可以作为指导流域科学评价的相互印证, 这些高风险区域也是未来土地规划的重点关注区域, 应加强生态风险管控^[30]。

多尺度生态风险评价有利于在整体生态风险评价的基础上, 发现高风险区域后进一步缩小研究区域范围, 精准识别风险区域。在本研究中, 流域尺度整体高风险区占比很小, 但在水源保护区尺度上, 特别是一级水源保护区, 高风险区域集中分布且占比较流域尺度相对较大, 体现了水源保护区的生态脆弱性, 也是未来流域治理和规划的重要区域, 应立足山水林田湖草系保护与治理, 围绕密云水库流域水资源保护、水生态修复, 构建以水源保护区水域空间为核心, 水陆交错带空间为缓冲屏障的水源保护格局, 推进流域水生态保护与可持续发展^[31]。以网格的形式对流域进行空间划分, 不仅可更好的表达空间异质性, 也有利于流域网格化精细化管理, 提高精准管控手段。但本研究仅局限于基于土地利用变化的相关研究, 并未考虑其他生态过程, 以流域为单元进行多尺度、多过程的生态模拟与优化调控是论文下一步需要探索的方向^[32—33]。

4 结论

本文以土地利用遥感数据为基础, 运用ArcGIS、Fragstats和GeoDa等软件对1990—2018年四个时期密云水库流域土地利用变化、景观指数变化、流域尺度和水源保护区尺度开展生态风险及空间自相关情况进行分析评价, 主要得到以下结论：

(1) 密云水库流域主要土地利用类型为林地和草地, 随着国家各项退耕还林政策的落实, 耕地减少了121.03km²；随着城市建设的扩张, 未利用地面积减少6.36 km², 建设用地面积增加200.56km², 2018年建设用地的面积是1990年建设用地面积的3.5倍。

(2) 景观格局指数显示, 景观斑块数(NP)、斑块密度(PD)、景观形状指数(LSI)、分离度指数(SPLIT)、香农多样性指数(SHDI)呈逐渐上升趋势, 蔓延度指数(CONTAG)、最大斑块指数(LPI)呈逐渐下降趋势, 景观指数变化体现了流域景观趋于局趋于复杂和分散, 破碎化程度加剧。

(3) 生态风险评价结果显示, 在流域尺度, 高风险区域面积逐渐向低风险区域转移, 呈现了“边缘高、中间低”的分布规律, 高风险区域集中分布在流域边缘的南部密云区、兴隆县、赤城县和流域北部丰宁县。在水源保护区尺度上, 呈现了“中间高、边缘低”的分布规律, 高风险区分布较为集中, 高风险区面积逐减少, 生态风险趋于减弱, 与全流域生态风险变化趋势相同。

(4) 空间关联指数结果显示, 流域景观生态风险呈正相关现象, Moran′s I指数均大于0, 1990年Moran′s I最大为0.322, 空间相关性最强, “热点”区域逐年增加, 局部空间自相关分布变化格局与景观生态风险分布变化趋于一致。