滨州市生态系统服务价值与生态风险时空演变及其关联性

李俊翰,高明秀

1 山东农业大学资源与环境学院, 泰安 271018 2 土肥资源高效利用国家工程实验室, 泰安 271018 3 湖南师范大学资源与环境科学学院, 长沙 410081

人口扩张与经济发展促使人类对地表空间占用范围和自然资源的开采力度不断增大,人类活动迅速而剧烈地改变着原始自然环境与生态系统,不可避免地造成大气、土壤和水资源的污染和破坏,反向影响人类自身福祉[1-3]。为了寻求高质量的生活品质和可持续的人类社会发展,人们愈发关注自身生活环境状况,人居环境、生态环境质量与安全的相关研究成为热点[4-5]。

生态系统服务评价和生态风险评价是生态环境评价的重要类型,与生态安全评价关系十分密切[6-7]。生态系统服务是人类从生态系统中获得的各种惠益,包括有形的物质产品供给与无形的服务提供两方面[8-9]。生态系统服务功能是实现生态系统安全的基本保证,从正面表征生态安全状况。生态风险是生态系统及其组分对外部干扰作出退化反映的可能性,指一定区域内具有不确定性的事故或灾害对生态系统的结构和功能可能产生的不利影响[10-11]。生态风险评价是对生态环境状况未来发展变化的预判,是风险控制的前提,从反面表征区域生态安全状况[12]。现阶段,对于二者的评价已由相对独立逐步走向融合,特别是基于生态系统服务的生态风险评价逐渐成为新的热点和研究方向[13-14]。曹祺文等[15]提出了基于生态系统服务的景观生态风险评价(ESRISK)框架,阐述了实施途径与优势。Dong等[16]也建立了类似的风险评估框架,基于甘孜地区土壤保持服务和防沙服务动态变化,评估生态风险及其变化,提高了对生态系统变化、生态系统服务和生态风险之间关系的认识。 刘世梁等[17]将生态系统服务功能的提升作为生态风险评价的重要内容,提出了基于土地整治全过程分析的生态风险识别、评价与管控思路。Liao等[18]将重要生态系统服务区作为风险受体,对厦门市城市综合生态风险进行综合评价和预测。宋科等[19]通过将生态系统服务价值的定量化与生态风险分析的数学模型相结合,建立了生态风险分析的EVR模型,定量化研究生态系统风险,并以浙江省舟山市为例验证了这种方法的可行性。由此可见,生态系统服务和生态风险评价的融合明显地表现为将生态系统服务引入到生态风险评价中,以提高生态风险评价的全面性和时效性[13,20]。

然而当前的研究中针对同一时期同一区域生态系统服务价值评估和生态风险评价的空间分布关系与相互作用机理研究较少。莫宏伟等[21-22]对农牧交错区榆林市和风沙过渡区神木县,宫继萍等[23]对黑河中游民乐县进行了土地生态系统服务价值评价和生态风险评价,但仅仅讨论了各自的空间分布和动态变化,缺乏对二者之间相互关系的深入探讨。鉴于此,该文以滨州市为研究区,核算2000—2015年生态系统服务价值和景观生态风险指数,从现状价值和潜在风险、正面和反面多角度描述区域生态环境现状及变化趋势,旨在基于生态系统服务价值和生态风险评价,探讨生态系统服务价值和生态风险指数在数量变化、空间分布和时序变化等方面的相互关联性,为构建区域生态安全格局,保护区域生态安全,实现可持续发展提供参考。

1 研究区概况

滨州市(36°41′—38°16′ N,117°15′—118°37′ E)位于山东省北部,黄河三角洲腹地,地处黄河三角洲高效生态经济区、山东半岛蓝色经济区和环渤海经济圈、济南省会城市群经济圈“两区两圈”叠加地带,辖区总面积9600km2,是黄河三角洲区域内最大的行政区。全境地势南高北低,主要分为低山丘陵和黄河冲积平原两大地形区,自然坡降比 1/4000—1/8000。滨州市属温带大陆性季风气候,夏热多雨,冬冷且长,四季分明,多年平均降水量575.4mm。受历史海陆变迁和地下水位高、矿化度大等因素影响,土壤盐渍化现象严重,境内盐碱地广泛分布。

2 材料与方法

2.1 数据来源与处理

选取2000(ETM+)、2005(TM)、2010(TM)和2015(OIL)年4—6月间的遥感影像数据(来自地理空间数据云平台http://www.gscloud.cn),利用ENVI5.2软件采用面向对象法进行土地利用分类,并基于2015年影像进行空间配准。参考全国遥感监测土地利用/覆盖分类体系和土地利用现状分类标准(GB/T 21010-2007),结合滨州市的自然环境与土地利用特征,划分为耕地、林地、建设用地、水域、盐田和未利用地等6种类型。分类后参考天地图和滨州市土地利用现状图等资料进行图像分类后处理,提高分类精度。每期分类结果随机选取1000个样点与实际调查地类进行精度验证,总体精度为82.82%、85.93%、86.97%和93.79%,Kappa系数为0.76、0.80、0.82和0.91,符合研究要求。

利用ArcGIS 10.2软件,采用3km×3km的单元网格对滨州市景观格局类型进行格网化重采样,最终形成1158个评价小区(图1)。基于评价小区计算生态系统服务价值和生态风险指数,并赋值于评价小区中心点,进行空间插值,分析空间格局分布特征。

图1 滨州市行政区及评价小区划分Fig.1 The evaluation cell on administrative area map of Binzhou city

2.2 生态系统服务价值核算

基于谢高地等[24]提出的生态系统服务价值化方法,定义1hm2面积农田生态系统粮食生产的净利润作为1个标准当量因子的生态系统服务价值量(标准当量)。本区域内耕地为旱地,主要种植小麦和玉米;林地以枣园、桑园等果园和针阔叶混交林为主;未利用地则为盐碱地和滩涂;盐田在区域东北沿海集中分布,多进行统一的生产、利用管理。本研究尺度较小,区域自然地理特征较为统一,农业耕作和养殖方法相同,区域各生态系统内部属性近似,可以应用当量赋值的方法进行服务价值核算。参考已有研究[25-26]并结合研究区实际状况进行标准当量的系数修订,修正模型如下：

Vi=(Swi×Pwi×ywi+Sci×Pci×yci+Sri×Pri×yri)×ei×fi

式中,V是研究期内滨州市标准当量(元/hm2);Vi是第i年的滨州市标准当量(元/hm2);n为研究期年份;Swi、Sci和Sri分别是第i年滨州市小麦、玉米和稻谷的播种面积占3种作物播种总面积的百分比;Pwi、Pci和Pri是第i年3种作物的全国平均净利润;ywi、yci和yri是第i年的单产修正系数;ei和fi是第i年的经济修和肥力修正系数。具体修正系数值计算公式如下：

fi=μfi×ρfi

参考曾杰等[27]的相关研究,单产修正值为滨州小麦、玉米和稻谷的平均单产值与当年全国单产值的比值。参考杨凤海等[25]的相关研究,αei是滨州人均地区生产总值与当年人均国内生产总值的比值,βei是研究期内历年滨州人均地区生产总值平均值与第i年滨州人均地区生产总值的比值。μfi是化肥修正系数,即滨州市平均单位面积农用化肥施用量(折纯量)与当年全国三种粮食作物平均施肥量(折纯量)的比值;ρfi是土壤生产力修正值,滨州市土壤生产力指数与全国平均土壤生产力指数的比值,参考孙振宁等[28]的研究结论,滨州市所处北方平原区土壤生产力指数为0.198,高于全国平均水平0.175,比值为1.13。

基础数据来自历年中国统计年鉴、山东统计年鉴、滨州统计年鉴和全国农产品成本收益资料汇编,相关数据经过基础数据计算得到。最终得到2000—2015年滨州市的标准价值当量因子平均值为1853.63元/hm2,作为研究期内滨州市标准当量基准值。

依据滨州市的实际情况,参考谢高地等[24]的研究成果,对各生态系统的服务价值系数进行修正：耕地对应旱地,林地对应阔叶林,建设用地为0,水域按水系和湿地1∶1的比例计算,未利用地按裸地和草甸1∶3的比例计算,盐田按建设用地、荒漠和冰川积雪3∶4∶3的比例计算。滨州市生态系统单位面积生态服务价值如表1所示。

表1 2000—2015年滨州市生态系统单位面积生态服务价值/(元/hm2)

2.3 景观生态风险指数核算

参考前人的研究成果[29-31],基于土地利用类型面积比重构建景观生态风险指数,定量描述评价小区的生态风险大小。

Ri=Ei×Fi

Ei=aCi+bNi+cDi

其中,ERIk为第k个样方的景观生态风险指数,Aki为样方内i类景观面积,Ak为k样方内景观总面积,Ri为i类景观损失度指数,Ei为干扰度指数,Fi为脆弱度指数,Ci为破碎度指数,Ni为分离度指数,Di为优势度指数。a、b、c为计算景观干扰度指数Ei的指标权重,参考前人研究[31-33],分别赋值0.5、0.3和0.2。景观脆弱度指景观类型抵抗外部干扰能力的大小,依据各类型的转出动态度排序,确定研究区景观类型的脆弱性排序,从高到低依次为：未利用地>林地>水域>耕地>盐田>建设用地,进行归一化处理后得到各景观类型的脆弱度指数Fi,分别为0.2857、0.2381、0.1905、0.1429、0.0952和0.0476。

2.4 灰色关联分析模型

灰色关联分析是依据序列曲线几何形状的相似程度来判断不同序列之间的联系是否紧密,灰色关联度是定量化描述各因素之间发展趋势相似或相异程度,适合动态的历程分析[34]。设序列X0,Xj长度相同,且初始值不为0,则两序列的灰色综合关联度为：

ρ0i=θε0i+(1-θ)r0iθ∈[0,1]

式中,ρ0i为两序列的灰色综合关联度,ε0i和r0i为两序列的灰色绝对关联度和灰色相对关联度。一般地,取值θ=0.5,θ取值越大越能体现绝对量之间的关系,θ取值越小越能体现变化速率之间的关系。相关指标应用灰色系统理论建模软件(GSTA V7.0)计算得到。

2.5 双变量空间自相关模型

空间自相关模型能够反映某种因素在空间位置的相关程度,分为全局空间自相关和局部空间自相关[35]。为了探究生态系统服务和生态风险的空间相关特征,采用双变量空间分析模型[36-38],利用全局自相关系数Moran′sI指数反映整体空间关联和差异状况,计算公式如下：

式中，Isr为单位面积生态系统服务价值s和生态风险指数r的双变量全局自相关系数;yi,s和yi,r为第i个评价小区单位面积生态系统服务价值和生态风险指数,σs和σr为方差。

为了全面具体地反映研究区各部分之间的空间关联性,采用LISA进行局部空间自相关分析,以表述局部集聚和离散效应。依据空间分布关系划分为高高聚集(H-H)、高低聚集(H-L)、低高聚集(L-H)和低低聚集(L-L)4个集聚类型区。

3 结果与分析

3.1 生态系统服务价值特征分析

2000—2015年滨州市生态系统服务价值核算结果和结构变化情况见表2。2000—2015年滨州市生态系统服务价值总量由168.06×108元持续增长到205.30×108元,年均增长率1.34%;单位面积生态系统服务价值由18553.00元/hm2增长到22663.34元/hm2。从研究区土地利用结构来看,各土地利用类型生态系统服务价值比例大小为水域﹥耕地﹥林地﹥未利用地﹥盐田﹥建设用地;耕地和未利用地生态系统服务价值持续减少,年均变化率分别为-0.93%和-7.04%,林地、盐田和水域生态系统服务价值持续增加,年均变化率为6.49%、4.81%和2.06%,并且所有土地利用类型生态系统服务价值增减速率逐渐降低,趋于平缓。从生态系统服务价值结构功能来看,单项生态服务价值比例大小依次为调节服务﹥支持服务﹥供给服务﹥文化服务;研究期间,各单项服务价值呈现持续增加趋势,年均变化率大小依次是调节服务(1.69%)﹥文化服务(1.61%)﹥支持服务(0.57%)﹥供给服务(0.38%),但是增长速率逐渐降低。

表2 2000—2015年滨州市生态系统服务价值结构变化表

通过计算评价小区单位面积生态系统服务价值,赋值于评价小区中心点,基于ArcGIS平台普通克里格方法进行空间插值,采用几何间隔分类,结合研究区实际状况,划分为低[0,8500)、较低[8500,10000)、中[10000,17000)、较高[17000,45000)、高[45000, ∞)共5个等级,形成4期单位面积生态系统服务价值空间格局分布图(图2)。滨州市单位面积生态系统服务价值主要以中和较高两个等级为主,面积平均占比为37.58%和35.54%;中和较低等级面积比例逐渐减少,高、较高和低级面积比例逐渐升高,但所有等级的面积变化速率都逐渐降低;北部平原、西南部丘陵岗地和南部黄河两岸区域单位面积生态系统服务价值较高,主要为高级和较高等级分布区域,中部为中等级分布区,西部和南部主要是较低和低级分布区域。

图2 2000—2015年滨州市单位面积生态系统服务价值空间格局Fig.2 The spatial pattern of per unit area ecosystem service values in Binzhou from 2000 to 2015

3.2 生态风险特征分析

滨州市生态风险评价结果见图3和图4。研究期间滨州市生态风险指数呈现缓慢降低的变化趋势。2000年、2005年、2010年和2015年滨州市生态风险指数值分别为0.0222、0.0215、0.0215和0.0209,总体下降5.86%,整体生态安全状况有所改善。从土地利用类型来看,耕地和未利用地生态风险指数值逐渐降低,研究期间共减少14.67%和33.91%,促使研究区生态安全状况得到一定程度改善。建设用地、林地、水域和盐田生态风险指数值呈现逐渐增加趋势,但增加速率逐渐降低,其中水域生态风险指数值增加最多,为46.40%,林地增加最少,为26.88%。

通过计算评价小区生态风险指数,采用几何间隔分类,结合研究区实际状况,划分为无[0,0.0120)、弱[0.0120,0.0180)、低[0.0180,0.0220)、中[0.0220,0.0370)、高[0.0370, ∞)共5个生态风险等级。从滨州市4期生态风险指数空间插值(图3)可以看出,滨州市生态风险等级呈现出总体连片、局部集聚的分布特征。研究期内主要以中生态风险和弱生态风险为主,期末面积占比达到34.90%和35.92%;高生态风险和无生态风险分布面积较小,期末占比仅为2.00%和4.53%。无生态风险区在研究区东北部集中分布,并且在此区域逐渐扩张;弱生态风险区在研究区中部偏西、西南部和东南部连片分布,少部分分布在无生态风险区外围;低生态风险区主要分布在滨州中部,并呈现面积逐渐减少的变化趋势,主要向中生态风险等级转化;中生态风险区主要分布在滨州市北部和黄河沿岸附近区域,但北部分布区面积逐渐扩大,黄河附近分布区域面积逐渐减小,向低生态风险等级转变;高生态风险区在研究期初主要分布于东北部沿海地区和西南部长白山山区,东北分布区逐渐向无生态风险等级转化并逐渐消失,西南部分布区面积逐渐减小,但依然是高生态风险区,同时滨州市中部逐渐形成高生态风险集聚区。

图3 2000—2015年滨州市生态风险空间分布格局Fig.3 The spatial pattern of ecological risk in Binzhou city from 2000 to 2015

如图4所示,滨州市各生态风险等级转入和转出变化总量均呈现逐渐减少趋势,总体生态风险等级变化状况逐渐趋于稳定状态。研究期间,低生态风险等级转出面积最多,共转出366790.91hm2,占转出总面积的42.78%;中生态风险等级转入面积最多,共转入274953.98hm2,占转入总面积的32.07%;无生态风险是转出和转入面积最小的生态风险等级,共转出7181.92hm2,转入40773.47hm2,占转出、转入总面积的0.84%和4.76%。从转化速率来看,转出面积平均变化率大小依次为高生态风险(65.73%)﹥低生态风险(44.88%)﹥中生态风险(27.77%)﹥弱生态风险(19.84%)﹥无生态风险(10.42%),转入面积平均变化率大小依次为高生态风险(53.23%)﹥无生态风险(44.49%)﹥低生态风险(37.14%)﹥中生态风险(31.24%)﹥弱生态风险(24.41%)。总体而言,研究期内滨州市生态风险等级动态变化表现为转入强度高于转出强度的特征,2005年以前等级变化较为活跃,等级间转化活动明显,2005年之后转化活动逐渐趋于缓和,并且变化强度逐步稳定。

图4 2000—2015年生态风险等级转入转出情况Fig.4 The situation of ecological risk grade transformation in Binzhou from 2000 to 2015

3.3 生态系统服务和生态风险的关联性分析

3.3.1数量相关关系

分别基于4个时期1158个评价小区的单位面积生态系统服务价值和生态风险指数值绘制散点图,并计算相关系数,进行函数拟合。如图5所示,数据点逐渐呈现出类似三角形的分布格局,单位面积生态系统服务价值和生态风险指数没有表现出明显的函数变化关系,仅在一定程度上表现为类似正相关关系,即生态风险指数值增大,单位面积生态系统服务价值也增大。基于SPSS软件进行Kolmogorov-Smirnov正态性检验,单位面积生态系统服务价值和生态风险指数两个变量都不服从正态分布,故不能采用Pearson相关分析,选用Spearmen相关分析两变量间的相关关系。4个时期Spearmen相关系数分别为0.743、0.772、0.736和0.757(P﹤0.01),表明单位面积生态系统服务价值和生态风险指数具有较强的正向等级相关性。选取线性、对数、倒数、二次、三次、复合、幂、S、指数和Logistic共10个函数模型进行函数关系拟合,计算拟合确定系数R2,得到2000、2005、2010和2015年R2值范围分别为0.156—0.389,0.183—0.414,0.174—0.383和0.083—0.420,函数拟合程度较差,说明二者没有明显的函数关系。

图5 2000—2015年单位面积生态系统服务价值-生态风险指数散点图Fig.5 The Scatter plot of unit area ecosystem service values and ecological risk index from 2000 to 2015

基于网格分析为基础,选取评价小区单位面积生态系统服务价值、单位面积供给服务价值、单位面积调节服务价值、单位面积支持服务价值和单位面积文化服务价值等5个指标与各样方生态风险指数值进行灰色关联分析,选取θ=0.5为基准值,计算灰色关联度。随后扩大和缩小θ取值,分析绝对量和变化速率的关系,仅选取θ=0.3、0.5和0.7时的灰色综合关联度计算结果进行展示(表3)。在取值范围内,不同θ取值的灰色综合关联度变化趋势相同,表明在绝对量关系和变化速率关系中各生态系统服务功能与生态风险指数的关联性基本一致。2000年单位面积供给价值与生态风险指数值的综合关联度最低,表明滨州市生态系统供给功能变化对区域生态风险变化的关联性最小;单位面积文化价值在4个功能价值中与生态风险指数值的综合关联度最高,表明滨州市生态系统供给功能变化对区域生态风险变化的关联性最大;单位面积生态系统服务价值与生态风险指数值的综合关联度最高,表明生态系统服务功能的综合作用对区域生态风险的影响作用大于各单项功能的影响作用。2005年之后,单位面积供给价值与生态风险指数值的综合关联度最高,高于单位面积生态系统服务价值与生态风险指数值的综合关联度,表明供给功能变化相比于生态系统综合服务功能变化更能体现生态风险的变化状况;单位面积文化价值与生态风险指数值的综合关联度最低,文化功能对生态风险的影响保持低水平状态。

表3 2000—2015年滨州市生态系统服务功能与生态风险指数灰色综合关联度

Table 3 The synthetic degree of grey incidence between unit area ecosystem service values and ecological risk index in Binzhou from 2000 to 2015

服务类型价值Value of service categoryθ=0.3θ=0.5θ=0.72000年2005年2010年2015年2000年2005年2010年2015年2000年2005年2010年2015年单位面积服务价值Unit area values of ecosystem services 0.83630.73820.64470.66270.74020.67010.60340.61620.64410.60210.56200.5697单位面积供给价值Unit area values of provisioning services0.66430.80010.72860.76940.61800.71450.66340.69250.57180.62880.59810.6156单位面积调节价值Unit area values of regulating services0.78570.71820.63250.65450.70410.65590.59460.61040.62240.59350.55680.5662单位面积支持价值Unit area values of support services0.79300.74660.66110.65270.70930.67620.61510.60910.62570.60580.56910.5655单位面积文化价值Unit area values of cultural services0.75020.69260.62390.63430.67890.63770.58860.59600.60750.58250.55330.5577

3.3.2双变量空间自相关

图6 2000—2015年单位面积生态系统服务价值与生态风险指数的全局Moran′s I指数(P<0.05)Fig.6 The global Moran′s I of unit area ecosystem service values and ecological risk index from 2000 to 2015(P<0.05)

采用双变量空间自相关模型,分析滨州市生态系统服务(单位面积生态系统服务价值)和生态风险(生态风险指数)的空间关联特征。结果表明(图6),2000年、2005年、2010年和2015年全局Moran′sI指数(p=0.05)分别为0.283、0305、0.269和0.285,表明生态系统服务与生态风险之间存在一定的正向空间关联性。从各分类服务价值来看,2000—2015年,支持服务价值和文化服务价值与生态风险指数的空间相关性高于综合服务价值与生态风险指数的空间相关性,供给服务价值和调节服务价值与生态风险指数的空间相关性则低于综合服务价值与生态风险指数的空间相关性;供给服务价值、文化服务价值与生态风险指数的全局Moran′sI指数变化趋势与综合服务价值的变化趋势相一致,调节服务价值、支持服务价值与生态风险指数的全局Moran′sI指数变化趋势一致;2000年、2005年和2010年支持服务成为区域主导服务,服务价值与生态风险的空间相关性主要反映其价值与生态风险的空间相关,2015年文化服务成为区域主导服务,空间相关主要反映区域文化服务价值与生态风险指数的空间相关性。

局部空间自相关LISA分析结果如图7所示,生态系统服务与生态风险局部关联性显著区域主要分布在北部无棣县、沾化区和南部邹平市,分为高价值-高风险(HH)、高价值-低风险(HL)、低价值-高风险(LH)、低价值-低风险(LL)和不显著(NS)5种模式。

图7 2000—2015年单位面积生态系统服务价值与生态风险指数空间自相关Fig.7 The Spatial autocorrelation between unit area ecosystem service values and ecological risk index in Binzhou from 2000 to 2015

高价值-高风险区：指单位面积生态系统服务价值高,并且生态风险等级高的区域。主要分布在滨州市西北部无棣县、东北部沾化区和西南部邹平县。北部集聚区主要集中于水域和未利用地交接处,水域能提供高价值的生态系统服务,却面临面积不断减小、向未利用地转化的高生态风险;西南集聚区以山区林地为主,是滨州南部生态系统服务供给中心,却因地表相对破碎、坡度较大,极易发生滑坡、泥石流等自然灾害。随着无棣县和沾化区推广冬枣种植,县域内枣林面积逐步增加,促使北部集聚区逐渐向南迁移,面积不断增加;南部则因山地开发,修建梯田,降低了生态系统服务价值和生态风险,集聚区面积逐渐减小并趋于稳定。

高价值-低风险区：指单位面积生态系统服务价值较高,但是生态风险较低的区域。研究期初主要分布在低价值-低风险区周边的水域和林地区域,后逐渐集中于北部盐田周边的水域分布区。相对于盐田,水域拥有更高的生态系统服务供给能力,但受气候的影响作用十分明显,极易向未利用地转化,是最有可能发生生态风险快速增高的区域。

低价值-高风险区：指单位面积生态系统服务价值较低,但是生态风险较高的区域。主要分布在高价值-高风险区的周边区域,即水域和林地周边的耕地区域。受高价值地类集聚效应及自身所处地类性质影响,所能提供的生态系统服务价值相对较低;受区域相对干旱的自然条件,在土壤盐渍化作用影响下,极易向荒地等未利用地转变,潜在生态风险较高。

低价值-低风险区：指单位面积生态系统服务价值低,生态风险低的区域。研究期初主要分布在无棣、阳信、邹平、博兴和滨城,北部集聚区以未利用地为主,其余集聚区以连片分布的耕地为主。研究期间,北部集聚区面积逐渐扩大,集聚地类向盐田转变,并趋于稳定;其余部分集聚区面积快速减小,并向建设用地集中。此类型集聚区主要分布在人为因素干扰强度大的区域,所提供的生态系统服务主要以人类最直接的需求为主,人类活动对区域生态环境的改造,使其处于相对长期稳定状态,生态风险较低。

不显著区：指单位面积生态系统服务价值与生态风险指数集聚关系不显著的区域。占比研究区面积的80.71%—84.10%,主要地类为成片分布的耕地和建设用地,生态系统服务价值和生态风险指数分布都较为平均,不存在明显的高值中心和低值中心。

3.3.3时间相关关系

基于1158个评价小区单位面积生态系统服务价值和生态风险指数值,将各评价小区2000年生态风险指数值与2000—2005年单位面积生态系统服务价值变化结果叠加,2005年生态风险指数值与2005—2010年单位面积生态系统服务价值变化结果叠加,2010年生态风险指数值与2010—2015年单位面积生态系统服务价值变化结果叠加,统计评价小区数目,探究各时段期初生态风险强度对期内生态系统服务价值变化的影响,结果如表4所示。

由表4可知,期初为高、中和低生态风险等级的评价小区,在3个研究时间段内单位面积生态系统服务价值表现为减少的数目都要高于表现为增加的数目,3个时段总计的生态系统服务价值增加数目占比分别为51.20%、50.57%和51.55%,都达到对应生态风险等级评价小区数量的一半以上;弱生态风险等级的评价小区单位面积生态系统服务价值变化主要表现为增加,占所有弱生态风险小区的47.01%;无生态风险等级的评级小区单位面积生态系统服务价值变化主要表现为不变和增加,分别占比41.48%和38.51%。总体来说,高等级生态风险区域预示着未来该区域生态系统服务价值最可能呈现减小趋势,低等级生态风险区域在未来生态系统服务价值往往表现为保持稳定或增加的变化趋势。

表4 2000—2015年滨州市生态风险强度对生态系统服务价值变化的影响

4 结论与讨论

(1)2000—2015年滨州市生态系统服务价值总量保持持续增长,由168.06×108元持续增长到205.30×108元,但增长速率逐渐降低并趋于平缓,年平均增长率1.34%;生态风险指数值持续降低,研究期间总共下降5.86%,总体生态风险等级逐渐趋于稳定;区域生态安全状况总体有所改善。

(2)在格网尺度,滨州市单位面积生态系统服务价值和生态风险指数具有较强的正向等级相关性,但没有明显的函数关系;单位面积供给服务价值与生态风险指数灰色综合关联度最高,同步变化程度较高,表明供给服务变化对该区域生态风险指数变化的影响作用较大。

(3)滨州市单位面积生态系统服务价值与生态风险指数之间存在一定的正向空间关联性。 2000—2015年,高价值-高风险区和高价值-低风险区面积不断增加,低价值-高风险区和低价值-低风险区面积逐渐减少;生态风险指数对生态系统服务价值变化具有一定的预示作用,高等级生态风险区域的生态系统服务价值存在极大的减小变化的可能性。邹平、无棣和沾化是高价值-高风险区的主要集聚区,成为全市生态系统服务的主要供给区和生态风险高危区。生态系统服务价值的不断提高表明生态环境的不断改善,较高的生态风险指数则揭示了生态系统结构和功能损伤的可能性,因此要特别注重此类区域的环境保护和生态安全建设的质量。

生态系统服务价值评价和生态风险评价的综合运用,能更好地将人类福祉与生态环境变化相联系,更好地反映区域社会发展所引起的生态环境改变情况,评价区域可持续发展现状及潜力,为区域生态环境保护和风险管理提供决策支持。由于该文进行的生态系统服务价值和生态风险指数核算都是基于区域土地利用类型(景观类型)结构所进行的,是两个结果之间的相互作用分析,尚不能很好地描述和反映生态系统服务与生态风险评价直接的内在联系;对二者相互关系的确定更多地是从相关性角度出发,对区域主导服务的识别和各单项服务与生态风险的权衡、协调研究只进行了初步尝试,在研究的广度和深度上还有较大空间。探讨生态系统服务与生态风险的内在机理,特别是利用不同的方法对二者进行服务价值和风险指数核算,从而进行有区分度的数量、空间相关性分析;构建生态系统服务-生态风险评价框架,提高二者综合评价对生态安全保护的有效性和时效性;理清生态系统服务和生态风险相互影响、作用和变化关系,特别是明辨单一生态系统服务对区域生态服务和风险的作用效果,应当作为下一步研究的重点和方向。