基于生态系统服务的城市化区域生态风险表征方法研究

刘长峰,侯 鹰,陈卫平,崔昊天

1 中国科学院生态环境研究中心城市与区域生态国家重点实验室, 北京 100085 2 中国科学院大学, 北京 100049

中国在过去40年经历了快速城市化的过程,在推动各个地区社会经济快速发展的同时,也对区域生态环境产生了诸多不利影响[1- 3],如人类生产生活产生的大量污染物损害了生态系统的结构、功能和健康[4],城镇建设用地扩张导致城市周边的森林、草地、湿地、农田等自然或半自然生态系统面积下降,景观斑块破碎,并造成生态系统服务损失[5- 6]。城市化导致的生态系统的不利变化目前已成为学界关注的热点问题[7]。

城市化的生态风险是指城市经济发展、城区建设和城市人类生产生活活动所导致的生态系统等的不利变化或对人类健康造成的不良影响的程度和可能性[8]。过去对于城市化的生态风险研究主要集中在土地利用变化和环境污染两个方面:(1)以区域景观格局作为生态风险评价终点,基于景观指数和土地利用类型得分构建生态风险指数,用以评估土地利用变化带来的生态系统受损的风险[9- 12]；(2)城市化过程中大量污染物的排放影响人体健康和区域生态系统特定功能的风险[13- 15]。以上研究的风险评价终点往往只是生态系统的某一种属性或过程,不能综合地反映城市化导致的区域整体的生态风险水平。

生态系统服务能够反映生态系统的整体性,使评价者能够明确保护对象和属性,其货币形式的评价结果便于同决策者交流,是理想的生态风险评价终点[16-17]。近年来以生态系统服务作为生态风险评价终点的研究逐渐增多,Thomsen 等[18]基于土壤生态系统服务构建了土壤污染的生态风险评估概念模型；Munns等[19]提出了以生态系统服务的损失耦合场地生态风险评估和自然资源破坏评估的概念框架；曹祺文等[20]提出了基于生态系统服务的景观生态风险评价框架。实例研究方面,Kang等[21]耦合区域生态系统服务价值和生态系统健康评估了京津冀城市群的生态风险；Pártl等[22]沿用传统风险定义从生态系统服务的供给脆弱度和风险源的胁迫程度两个方面,评估了捷克共和国的综合生态风险水平。然而,目前基于生态系统服务的生态风险评价研究仍处于起步阶段,有关风险表征方法的研究还很缺乏。此外,由于缺乏同生态系统服务相关的风险阈值,使得传统的风险商表征方法不能适用。因此,大多数现有研究还停留在评价框架或概念模型阶段,实例研究较少。

针对已有研究的不足,本研究旨在建立适用于城市化区域的基于生态系统服务的生态风险表征方法。将多种类型的生态系统服务纳入到生态风险评估中,参考TOPSIS (Technique for Order Preference by Similarity to an Ideal Solution) 模型(优劣解距离法),以空间上各类生态系统服务价值的最大损失量作为风险参考值,构建了生态风险综合指数,建立了以生态系统服务为评价终点的生态风险表征方法,并使用构建的方法分析了北京市的生态风险；最后使用基于蒙特卡洛模拟的不确定性分析和参数敏感性分析对生态风险指数的可靠性进行了分析。

1 研究方法

1.1 生态风险评估框架

生态风险评估框架主要包括问题形成、风险分析和风险表征三个过程。本研究关注的是城市化区域生态系统服务下降的风险,城镇土地扩张是该类区域中城市化过程最显著的特征。过去很多研究表明土地利用变化是导致生态系统服务变化的最主要因素[6,23],因此城市化区域的风险胁迫主要是城市化带来的土地利用等人类活动,相应的风险评价终点是区域生态系统服务,服务下降是风险胁迫产生的不利生态效应(图1)。

图1 基于生态系统服务的城市化区域生态风险评价基本框架Fig.1 A basic ecological risk assessment framework for urbanized areas based on ecosystem services

1.2 生态风险表征方法

传统的评价环境污染物导致的生态风险是通过量化风险胁迫的程度(如污染物浓度)或危害程度(如健康损害指数)来表征生态风险水平[13- 15]。本研究中风险胁迫(土地利用)与风险评价终点(生态系统服务)之间的关系明确,因此以生态系统服务损失量来表征生态风险(图1)。由于缺乏关于生态系统服务损失的风险阈值,本研究提出以生态系统服务价值量的相对损失量来表征生态风险,构建适用于城市化区域的生态风险综合指数。

生态风险综合指数参考了TOPSIS模型[24]。首先找出区域内所有空间样本(栅格像元或矢量多边形区域等)中各类生态系统服务价值达到最大值和最小值的两个样本作为理想样本点(正理想样本点和负理想样本点),以正理想样本点的服务价值量作为风险阈值,以这两个理想样本点的欧氏距离作为该研究区域生态系统服务损失达到最大的表征,则综合风险指数定义为每一个真实样本与正理想样本点间的欧氏距离与生态系统服务损失达到最大的欧氏距离之比。具体计算方法如下:

(1)基于现有的某种生态系统服务价值计算方法计算出m个空间样本(栅格像元或矢量多边形区域等)的n类生态系统服务价值,提取所有空间样本构建生态系统服务价值矩阵ESV=(esvij)mn。

(1)

(2)为了消除量纲和数量级的影响,对生态系统服务价值矩阵ESV进行归一化处理(式2),之后计算加权样本矩阵Y=(yij)mn(式3)。

(2)

式中,rij:样本i的第j种生态系统服务价值的归一化数值；max(esv·j):第j种生态系统服务价值所有样本中的最大值；min(esv·j):第j种生态系统服务价值所有样本中的最小值。

(3)

式中,yij:样本i的第j种生态系统服务价值的加权归一化数值；ωj:第j种生态系统服务价值的权重,对不同生态系统服务类型的赋权体现了在生态风险评价中对不同生态系统服务的侧重不同。

(3)正理想样本点和负理想样本点的识别

(4)

(5)

式中,y+:正理想样本,即所有的生态系统服务价值量均达到最大值的理想样本点；y-:负理想样本,即所有的生态系统服务价值量均达到最小值的理想样本点；max(y·j):生态系统服务j在所有区域中加权归一化最大值；min(y·j):生态系统服务j在所有区域中加权归一化最小值。

(4)生态系统服务价值相对损失量的表征:

(6)

(7)

(5)风险指数的计算:

(8)

式中,Ri:样本i的生态风险指数,数值介于0到1之间。越接近于0表明i的生态系统服务相对损失量越小,即生态风险越低；越接近于1表明i的生态系统服务相对损失量越大,即生态风险最高。

(6)风险等级划分

为了减少风险等级划分的主观性,通过对样本性质进行探索,找到相应的风险等级分界点。由于构建的风险指数只有单一属性列表,故采用高斯混合模型对数据的聚类数目进行探索[25]。如果所有风险指数样本呈多峰分布且每一个峰均服从正态分布,则可以使用每一个正态分布峰的均值或分位点作为风险等级的分界线；如果不是正态分布,则使用Jenks自然断点分类法划分生态风险等级[26]。

2 案例研究

2.1 研究区概况

近几十年来北京市的城市化发展迅速,城市区域快速扩张,现已形成五条环城公路。一般认为五环以内的区域为北京的市区,五、六环之间的区域为北京的郊区；六环路以外的区域为北京的农村地区。北京市的人工地表主要集中于六环以内,人工地表率随着离城市中心的距离增加而明显下降(图2)。因此,北京市是典型的快速城市化区域,是研究快速城市化发展导致的区域生态风险的理想区域。

图2 2015年北京市土地利用类型分布图Fig.2 The land use map of Beijing in 2015

2.2 生态系统服务价值核算方法及权重的确定

生态系统服务价值核算方法主要有基于生态系统服务功能量的价值定量方法[27- 29]和当量因子法[30-34]。本研究使用当量因子法进行生态系统服务价值核算,当量因子法所需参数和数据量较少,易于在区域尺度核算各类生态系统服务价值[23,35-36]。基于当量因子法的生态系统服务价值量计算如下:

(9)

式中,E:1 个标准当量因子的生态系统服务价值量 (元/hm2),指 1 hm2全国平均产量的农田每年自然粮食产量的经济价值；Aik:区域i内土地利用类型k的总面积(hm2),本研究用到6种土地利用类型:森林、草地、耕地、湿地、人工地表和裸地；Djk:土地利用类型k所提供的生态系统服务j的当量因子,本研究参考了谢高地等制定的适用于全国陆地生态系统服务的当量因子表[34](表1),根据北京市的生态系统服务当量因子是全国平均水平的1.04倍,对应用于北京市的当量因子进行修正[31]。

本研究选取了适用于北京市的9类生态系统服务,使用熵权法对服务的权重进行计算[37- 38](表1)。

表1 生态系统服务单位面积当量因子和权重

2.3 不确定分析与参数敏感分析

由于本研究提出的生态风险指数无法使用实测数据验证,因此使用不确定性分析和参数敏感性分析对风险指数的可靠性进行评价。如果模型(本研究中为生态风险指数)计算结果的不确定性太大,说明模型的可靠性低[39- 40]；反之,如果模型模拟结果的不确定性可接受,则说明模型较为可靠。模型参数的敏感性分析是用于评价参数对模型输出结果影响程度的方法[41- 42],通过参数敏感性分析,可以找到对模型计算结果影响较大的参数,从而有针对性地对模型进行优化。

为提高风险指数可靠性评价的效率,不确定分析与参数敏感分析在模拟栅格中进行。模拟栅格的土地利用情况为森林、草地、耕地、湿地和人工地表各占20%,即不考虑栅格内部土地利用类型比例的不确定性,针对9种生态系统服务的权重和当量因子(当量因子的不确定性综合考虑5类土地利用类型当量因子之和)进行不确定分析与参数敏感分析。

2.3.1不确定性分析

本研究使用蒙特卡罗模拟的方法进行不确定性分析。蒙特卡罗模拟是通过从模型参数的概率分布中随机抽样,计算出一系列运算结果样本的过程[43]。本研究使用蒙特卡洛模拟结果样本的变异系数、最大变异范围和相对变异指数来表征风险指数计算结果的不确定性,三个指标值越小,说明风险指数计算结果的不确定性越小。

(1)变异系数(coefficient of variance,CV)

(10)

式中,μ为样本均值；s为样本标准差。变异系数是衡量一系列数据变异性的常用指标,适用于任何概率分布。

(2)最大变异范围

(3)相对变异指数(relative variability index,RVI)

2.3.2参数敏感性分析

本研究采用局部参数敏感性分析,检验每个参数值增加50%对风险指数计算结果的影响[42],计算每个参数的敏感性在全部参数的敏感性中的贡献度,计算公式如下:

(11)

式中,S(pi)为参数pi的敏感性贡献度,以某个参数变化导致风险指数计算结果的变化量绝对值占所有参数变化导致风险指数计算结果的变化量绝对值的总和的比例来表示；pi为参数i的实际取值；f(·)为模型运算公式。

2.4 数据与分析

本研究主要使用的数据为:2015年北京市土地利用分类30 m分辨率栅格数据[44],北京市行政边界数据(http://www.resdc.cn/Default.aspx),生态系统服务价值当量因子表[34]。 在ArcMap 10.5中将北京市划分为1 km×1 km的网格,然后使用局域统计工具提取出每个网格的土地利用配置；将数据导出后,在MATLAB 2017b中执行生态风险指数的运算以及蒙特卡罗模拟；使用R 3.4.3的mixtools工具包进行高斯混合模型训练以及数据的分布检验。使用ArcMap 10.5将生态风险计算结果进行空间化表达；使用MATLAB 2017b绘制风险指数频率分布直方图和蒙特卡罗模拟结果图；使用Excel 2016绘制参数敏感性分析结果图。

3 结果与讨论

3.1 生态风险等级划分结果

2015年北京市生态风险指数空间样本的频率分布直方图显示风险指数呈多峰分布(图3)。经过高斯混合模型的训练,将所有样本分为4个峰的贝叶斯信息准则值(Bayesian Information Criterion, BIC)和p值(Kolmogorov-Smirnov test)最高(表2),表明最合适的划分是分成4个子分布。进一步对每个子分布进行正态检验(Shapiro-Wilks test),结果表明4个子分布均服从正态分布(表3)。

图3 生态风险指数的频率分布直方图以及正态分布拟合曲线Fig.3 The frequency histogram and normal distribution fitting curve of ecological risk index

表2 高斯混合模型训练结果

表3 子分布正态性检验结果

生态风险指数的多峰分布反映了空间样本的分类,每类样本具有一定的相似性。基于这种分类可以在没有现成标准的情况下,比较客观地划分生态风险等级。本研究以4个子正态分布的均值为划分界限将2015年北京市的生态风险划分为五个等级(表4)。由于本研究提出的生态风险指数反映的是区域的相对生态风险,因此,风险等级也反映了风险的相对高低,当研究区域发生变化时,风险等级划分界限也会相应变化。

表4 生态风险等级划分以及相应风险指数范围

3.2 生态风险的统计描述及空间分布

2015年北京市生态风险指数最小值为0.19,最大值为0.99,平均值为0.64,生态风险总体上处于低风险接近中等风险水平。具体而言,极低生态风险区域面积为869.7 km2,仅占北京市总面积的5.3%,主要分布于北京市北部的部分区域。 低风险区域占据了最大的面积,为8041 km2,占北京市总面积的49%,主要分布于北京的西部和北部；中等风险区域和高风险区域面积分别为3889 km2和1640.9 km2,分别占北京市总面积的24.7%和10%,主要分布于北京市的东部和南部的郊区地带；极高风险水平区域面积为1641.1 km2,占北京市总面积的10%,主要分布于北京中部的城市区域(图4)。从整体上看,2015年北京市生态风险呈东高西低和南高北低的空间分布特征,并且具有明显的空间异质性和城郊梯度特征(图4)。已有其他区域的研究发现了相似的现象,如张小飞等[9]和杨勇等[11]的研究分别发现淮北城市综合生态风险和西安长安区土地利用生态风险具有很高的空间异质性；汪翡翠等[10]的研究显示京津冀区域中非城市地区的土地利用生态风险低于城市地区。此外,本研究发现北京市生态风险从高风险到低风险的空间变化趋势与该区域人工地表从密集到稀疏的空间变化趋势大体上一致(图2和图4),这表明城市化过程驱动的城镇扩张会导致周边区域生态系统服务下降,造成生态风险水平跨级增长,该结果与汪翡翠等[10]和Kang等[21]的研究结果相似。这两项研究表明京津冀的快速城市化地区较高和高生态风险等级区域范围明显增加。本研究案例的风险评价结果与过去多项研究结果的相似性表明本研究提出的风险表征方法合理可靠。

图4 2015年北京市不同生态风险等级的空间分布 Fig.4 Spatial distribution of different ecological risk levels in Beijing in 2015

3.3 生态风险指数与生态系统服务价值间的定量关系

为分析本研究提出的生态风险指数能否准确地反映区域生态系统服务的损失,以研究区域全部栅格为样本,进行生态风险指数与生态系统服务价值当量因子总数间的回归分析,结果显示生态风险指数的对数值与生态系统服务当量因子总数之间呈极显著线性关系(P<0.001,R2=0.9715)。因此,本研究提出的风险指数增加量可以反映生态系统服务价值的损失量。例如,当风险指数从0.4上升到0.5时,相应的生态系统服务价值当量因子数损失量为562个/km2。参考谢高地等[45]的研究,相应的生态系统服务价值损失量为191.44万元/km2。

3.4 不确定性分析和参数敏感性分析结果

3.4.1不确定性分析结果

本研究对生态风险指数进行了10000次蒙特卡罗模拟,正态性检验表明模拟结果呈正态分布(Shapiro-Wilks test,P>0.05)。模拟结果的变异系数仅为0.079,且3倍的标准差仅为0.081,相对变异指数(RVI)仅为-5%—5%。几个表征不确定性的统计指标值均较小,表明本研究提出的生态风险指数的不确定性是可以接受的,能够较好地表征一个土地利用类型组成确定的研究区域的生态风险(表5和图6)。

图5 生态风险指数与生态系统服务价值间的定量关系Fig.5 Quantitative relationship between ecological risk index and ecosystem service value

图6 蒙特卡罗模拟结果的累计概率分布 Fig.6 Cumulative probability distribution of Monte Carlo simulation results

表5 蒙特卡罗模拟结果的统计描述

3.4.2参数敏感性分析结果

基于蒙特卡罗模拟的结果,对每个参数的敏感性进行了分析,结果如图7所示。敏感性贡献度最大的三类参数分别是水文调节服务当量因子、文化娱乐服务当量因子和气体调节服务当量因子,贡献度分别为17.65%、13.06%和12.76%。9个权重的敏感性贡献度之和为17.33%,相比于所有当量因子,权重的敏感性贡献度较低。

图7 不同参数的敏感性贡献度比较Fig.7 Comparison of the contributions of sensitivity of different parameters食物权重:食物供给服务权重；原材料权重:原材料供给服务权重；气体权重:气体调节服务权重；气候权重:气候调节服务权重；水文权重:水文调节服务权重；环境权重:环境净化服务权重；土壤权重:土壤保持服务权重；生物权重:生物多样性维持服务权重；文化权重:文化娱乐服务权重；食物当量:食物供给服务当量因子；原材料当量:原材料生产服务当量因子；气体当量:气体调节服务当量因子；气候当量:气候调节服务当量因子；水文当量:水文调节服务当量因子；环境当量:环境净化服务当量因子；土壤当量:土壤保持服务当量因子；生物当量:生物多样性维持服务当量因子；文化当量:文化娱乐服务当量因子

图8 五类参数的敏感性贡献度比较 Fig.8 Comparison of the contributions of sensitivities of 5 types of parameters服务权重:各项生态系统服务权重；供给当量:供给服务当量因子；调节当量:调节服务当量因子；支持当量:支持服务当量因子；文化当量:文化服务当量因子

本研究中使用的生态系统服务价值定量方法为当量因子法,该方法所需参数少,但未考虑小尺度上生态过程的差异,相比于其他生态系统服务价值核算方法,在应用于中小尺度或人工生态系统中的结果不确定性更大[46]。本研究的不确定性分析结果表明即使使用不确定性较大的当量因子法,研究区域的生态风险评价结果的不确定性并不高。因此未来的研究可以进一步优化生态系统服务核算方法以获取更加准确可靠的风险评价结果。为了增加风险评价结果的准确性,需要考虑评价空间尺度、区域社会经济差异和自然环境差异等要素选取合适的生态系统服务价值计算模型。

4 结论

以多种类型生态系统服务作为城市化区域生态风险的评价终点,参考TOPSIS模型构建了基于服务价值量的生态风险表征方法,以北京市为例进行了方法的应用和不确定性分析。研究发现2015年北京市生态风险总体上处于低风险接近中等风险水平,生态风险具有明显的沿乡村到城市的梯度逐渐上升的特征,表明北京市城市区域的扩张造成周边区域生态系统服务的下降,导致生态风险水平的升高。提出的风险指数用以计算土地利用结构确定的区域内的生态风险时,计算结果变异较小,表明其具有较高的可靠性。此外,生态风险指数同生态系统服务当量因子间的定量关系明确,因此可用于快速估算生态系统服务价值损失量。为优化本研究提出的生态风险指数,未来的研究需要选择可靠性更高的生态系统服务价值量计算方法,以获得更加准确的风险评价结果。