基于生态服务价值的沱江流域土地景观生态风险时空分异研究

陈 丹， 蒋贵国， 张 妍， 涂 月， 邓元杰

(1.四川师范大学 西南土地资源评价与监测教育部重点实验室， 四川 成都 610068；2.四川师范大学 地理与资源科学学院， 四川 成都 610068; 3.西北农林科技大学 经济管理学院， 陕西 杨凌 712100)

生态风险是指一个种群、生态系统或整个景观的正常功能受外界胁迫，从而在目前和将来减少该系统内部某些要素或其本身的健康、生产力、遗传结构、经济价值和美学价值的可能性[1]。近年来，随着经济社会的快速发展，城市化发展迅速，使得人类活动更为复杂，产生了巨大的开发建设需求，这些都影响着土地利用景观格局，景观趋于破碎化，景观结构趋于复杂化，进而对土地利用的生态服务价值产生了一定的影响。复杂的自然过程和人类活动产生了一定的生态风险，威胁着人地关系的和谐[2]。区域生态风险评价是从区域尺度上描述和评价环境污染、人为活动或自然灾害对生态系统结构和功能等产生不利作用的可能性和危害程度[3]。近年来，对生态风险的研究已成为生态学、地理学等领域的热点话题[4-12]。

流域是一个是以河流为中心，集自然资源、社会、经济于一体的重要地理单元[13]，流域内蕴藏着丰富的自然资源和得天独厚的生态环境，孕育着众多的人口，承载着多种经济方式的发展，是土地利用类型中水域和其他土地利用类型交互作用的区域。不同土地利用方式的交互作用对流域的影响具有累加性，这种累加性达到一定程度会影响生态环境的结构和功能，进而影响生态系统提供服务的能力，产生一定生态风险。2015年5月，国务院发布《关于加快推进生态文明建设的意见》，其中流域生态环境质量的改善被列为目标之一。流域生态环境是否安全、是否遭受生态风险成为攸关民生的关键问题。因此，对流域进行生态风险评价并保护流域生态安全是很有必要的。

生态系统服务是指生态系统形成和所维持的人类赖以生存和发展的环境条件与效用，生态系统的结构、过程和功能直接或间接得到生命活动所需的产品和服务[14]。人类活动改变着土地利用景观格局，引起生态系统的结构和功能，进而影响生态系统提供服务的能力。不合理的人类活动会影响生态服务功能，引起生态风险，影响人类生产生活，最终阻碍可持续发展。

在已有的研究中，大多直接在土地利用数据基础上计算景观格局指数，由此计算生态风险指数，却少有基于生态服务价值的生态风险的研究。本文拟在沱江流域2000,2005,2010,2015年土地利用数据的基础上，计算生态服务价值，构建基于生态服务价值的沱江流域生态风险评价体系，明确土地利用景观生态风险评价的终点是保护生态服务功能，提升生态服务价值，并且将人类活动、景观格局、生态过程、生态服务功能联系起来。以期为基于生态服务价值的生态风险管控决策提供依据。

1 研究区概况

沱江流域(103.68°—105.76°E，28.87°—31.70°N)位于四川省中东部，包括了成都市、德阳市、资阳市、眉山市、内江市、自贡市、泸州市7个市级行政单位的25个县级行政单元(见附图4)，总面积26 210.56 km2。沱江全长712 km，地势总体上西北高东南低，境内有岷山、龙门山、龙泉山等山脉。据资料显示，沱江流域境内平均森林覆盖率仅为6.1%，是四川省的所有河流森林覆盖率的最低值。境内工业发达，如德阳的重型机器制造业、泸州的化工业等。截至2015年，流域境内总人口1.70×107人，GDP为7.01×1011元，人均GDP达4.12×104元/人，相比于2000年的4 800元/人，增涨了近8倍[15]，经济发展迅速。沱江流域近年来随着经济方式的多样化发展，城市化进程加快，城镇建成区快速扩张，资源紧缺、环境污染、生态脆弱等环境问题愈发明显，影响了生态系统提供服务的功能，生态环境遭受着一定的风险。沱江作为长江的重要支流，沱江流域作为长江经济带西部的重要组成部分，其生态系统遭受危险，不仅会威胁流域内人们的生产生活，还会制约下游的社会经济发展。

2 数据来源与方法

2.1 数据来源与处理

研究所用的数据及来源主要有：①行政边界矢量数据为全国1∶100万基础地理数据库，来源于国家基础信息中心(www.webmap.cn)。②DEM数据为ASTER-GDEMV2 30 m分辨率的数字高程数据，来源于NASA(search.earthdata.nasa.gov)。③归一化植被指数(NDVI)数据为中国500 m NDVI月合成产品(MODND1M)，来源于美国国家航空航天局(lpdaac.usgs.gov)。研究借助ENVI5.1平台，运用MVC最大值合成法，合成2000，2005，2010，2015年的NDVI年最大值。④净初级生产力(NPP)数据，为MOD17A3H，来源于美国国家航空航天局(lpdaac.usgs.gov)。⑤2000，2005，2010，2015年4期土地利用现状数据，分辨率为30 m，来源于中国科学院资源环境科学数据中心(http:∥www.resdc.cn)，重采样后得到分辨率为100 m的土地利用数据。⑥社会经济数据来源于四川省统计年鉴。参照土地利用现状分类标准(GB/T 21010-2017 )，将土地利用数据分为耕地、林地、草地、水域、建设用地、未利用地6类。

2.2 研究方法

2.2.1 风险小区的划分 采用网格法划分风险小区。综合考虑研究区域沱江流域的实际情况，研究以8 km×8 km的网格对研究区进行等间距采样，得到505个风险小区。计算各个网格中心的生态风险指数，以此代表该风险小区的生态风险水平。

2.2.2 基于生态服务价值的生态风险评价体系的构建 以往的研究多用景观指数和概率指数直接构建生态风险指数，但这种方法缺乏对特定生态过程的关注，且概率指数多由经验而得，主观性过大。生态系统服务是景观格局和生态过程相互作用的产物，将生态服务价值引入生态风险评价可以在较大程度上弥补上述缺憾[2]。本研究采用基于生态服务价值的生态风险评价方法，构建基于生态服务价值的沱江流域生态风险评价框架(见图1)。

图1 基于生态服务价值的沱江流域土地景观生态风险评价框架

(1) 损失度指数的计算。框架从损失度指数和概率指数两方面对生态风险进行综合表征，其中损失度由沱江流域土地利用各项生态服务价值的损失归一化处理得到，加权求和得到损失度，然后再用修正系数(公式5)对损失度进行修正，得到修正后的损失度指数。

①农田生态系统的经济价值计算方法。计算公式[16-17]如下：

(1)

式中：Pa——单位农田生态系统提供食物生产功能的经济价值(元/hm2)；i——作物种类；Pi——i种粮食作物的平均价(元/kg)；qi——i种粮食作物单产(kg/hm2)；mi——i种粮食作物面积(hm2)；M——粮食作物总面积(hm2)。

②沱江流域土地景观生态服务价值。研究采用以下公式[16,18]计算沱江流域土地景观生态服务价值：

VES=∑(Ak×VCk)

(2)

ESVf=∑(Ak×VCfk)

(3)

ESVa=ESV×R

(4)

式中：VES——生态服务价值； VCk——生态服务价值系数(元/hm2)；Ak——k种地类面积((hm2)； ESVf——单项生态服务价值(元)； VCfk——单项服务功能价值系数(元/hm2)；R——生态服务价值的修正系数。

研究采用修正系数R，来修正生态服务价值。VC的计算公式[19]如下：

(5)

式中：fgi——第i个格网的植被覆盖度，研究拟用归一化植被指数NDVI来表征植被覆盖度；fgmean——植被覆盖度的平均值； NPPi——第i个格网的净植被初级生产力； NPPmean——净植被初级生产力的平均值。

③基于格网的生态服务价值评价方法。本研究采用网格法(风险小区法)进行土地景观服务价值的评价，计算每个格网单元生态系统服务价值指数(ESVI)，作为格网中心点的生态服务价值水平，计算公式[20]如下：

(6)

式中：ESVIk——第k个格网生态服务价值指数(元)；Ak——格网单元总面积(hm2)；Aki——第k个格网中第i类土地类型的总面积(hm2)；Ei——第i类土地利用类型的生态服务价值系数(元/hm2)。在得到格网生态服务价值的基础上，计算各网格生态服务价值的差值，最后归一化处理得到每个网格的损失度指数L。

(2) 概率指数的计算。参考相关研究[2]，概率指数由概率评价指标体系(表1)得到，研究从地形、人为胁迫、生态恢复力、景观脆弱性4个维度，选取坡度、建设用地比例变化、净初级生产力、植被覆盖度、景观连通重要性、景观易损度6个因子，对沱江流域土地利用景观生态风险的概率进行量化，区别于仅由经验而得的概率值。

表1 基于生态服务价值的沱江流域生态风险评价的概率指标体系

对于景观格局脆弱度[21]，计算公式如下：

LVI=LSI×(1-LAI)

(7)

(8)

Ui=aFNi×bFDi×DOi

(9)

LAI=PRD×SHDI×SHEI

(10)

式中：LVI——景观格局脆弱度指数； LSI——景观敏感度指数； LAI——景观适应度指数；Ui——景观干扰度指数，表示外界因素对景观格局的干扰程度，由破碎度(FN)、分维数倒数(FD)、优势度(DO)加权求得，a，b，c——权重，分别取值0.5，0.3，0.2[22]。

(3) 基于生态服务价值的生态风险指数的计算。

计算公式如下：

ERI=Li×Pi

(11)

式中：ERI——生态风险指数；Li——第i个网格(风险小区)的损失度指数；Pi——第i个网格(风险小区)的概率指数[2]。考虑到生态风险值计算结果的聚集程度，将生态风险值划分为5个等级：低生态风险(ERI≤0.05)，较低生态风险(0.050.20)。

2.2.3 空间自相关分析法 研究采用空间自相关分析法分析生态风险的空间分布特征。

(1) 全局空间自相关的计算公式[23]如下：

(12)

(2) 局部空间自相关的计算公式[23]如下：

(13)

式中：xi，xj，wij和x的含义同公式(12)。表示相似值的空间聚集，LISA大于0表示该区域单元周围存在空间聚集现象。

2.2.4 半方差分析法 地统计软件GS+可以实现半变异函数理论模型的最优拟合，利用半方差的变异函数，从其中的主要参数可以定量分析生态风险在空间上的分布规律。计算公式[24]如下：

(14)

式中：r(h)——变异函数；h——样本空间距离(步长)；N(h)——样本对总数；Z(xi)，Z(xi+h)——空间位置xi，xi+h处的生态风险值。在理论拟合最优的基础上，利用普通克里格法编制不同时相沱江流域生态风险值的分布图。

3 结果与分析

3.1 土地利用动态变化分析

如表2所示，2000—2015年，沱江流域土地利用动态变化明显。其中，建设用地由548.94 km2增长到了1 077.75 km2，增幅达96.33%；林地面积减少473.60 km2，年均减少31.57 km2；耕地减少了150.20 km2；未利用地面积成倍数增长，但在土地利用总体变化中便显得并不明显。总体来说，沱江流域在2000—2015年间，建设用地面积扩张最多，土地利用转移以耕地转建设用地、林地转耕地为主。究其原因主要是，2000—2015年研究区经济快速发展，人均GDP由2000年的4 800元/人增加到了2015年的41 200元/人，年均增长率达50.48%。成都市作为中国西部的中心，发展迅猛，逐渐成为西南部甚至西部的经济增长极。经济的发展产生了巨大的开发建设需求，所以建设用地不断增加，耕地转变为建设用地，林地被开垦为耕地。

表2 2000-2015年沱江流域土地利用转移矩阵 km2

3.2 生态服务价值的变化分析

如表3所示，2000—2015年，沱江流域生态服务价值总量先增后降，耕地、林地对生态服务价值的贡献最大。除2000年外，2005—2015年耕地对当年生态服务价值的贡献率保持在35%以上，且呈上涨趋势，这是因为耕地承担的供给功能比较大，虽然耕地面积减少了，但是随着经济的发展，耕作技术的进步，粮食单产提高，耕地对当年生态服务价值的贡献率一直保持在较高水平。林地对生态服务价值的贡献率一直较高，但贡献率总体下降，由2000年的62.75%下降到了2015年的37.20%，这主要与林地面积减少有关，林地面积在研究期间减少了473.60 km2，面积变化率达11.82%。其中，建设用地的生态服务价值为负，计算过程如下：参考谢高地等[17]最新改进过的ESV计算方法及其制定的中国陆地生态系统单位面积生态服务价值当量表，并参考郑江坤等[25]得出的密云水库集水区的不同地类单位面积生态服务价值表和该地区的修正参数1.024，二者相除得到全国的建设用地单位面积ESV，再乘以四川省的修正参数1.35，得到四川省建设用地生态服务价值，将此作为最新的四川省建设用地生态服务价值指标，再运用公式(1)—(6)，得到沱江流域修正后的生态服务价值。

表3 2000-2015年沱江流域生态服务价值及变化 108元

3.3 生态风险的空间自相关分析

借助GeoDa软件平台，计算出沱江流域2005，2010，2015年的生态风险Moran’sI指数分别为0.714 669，0.706 358，0.736 719，并通过p=0.05的显著性检验。表明沱江流域生态风险指数呈显著的空间自相关关系，呈聚集分布，并非随机分布；且Moran’sI指数呈增加趋势，说明生态风险集聚现象越来越明显。“高高”、“低低”集聚的散点较多，说明生态风险高值与周围的高值、低值和周围的低值集中，同质性较为显著(见图2)。

由图2，可看出2005—2015年沱江流域生态风险值“冷点”(低—低)和“热点”(高—高)分布的具体位置。“热点”区域主要分布在流域上游的成都市、德阳市的大部分地区，包括龙泉驿区西部、彭州市北部、广汉市、什邡市北部、绵竹市北部。成都、德阳作为四川省的老牌经济发达城市，其经济地位在全省均处于前列，各经济主体活跃，经济发达且对发展经济的要求强烈，城镇建成区面积大且扩张迅速，建设用地比例高，生态服务价值低，生态风险高；就自然方面来说，龙门山、龙泉山横亘于成都、德阳所在地，地势起伏大，生态系统脆弱，易产生生态风险。

相比于“热点”区域，“冷点”区域分布相对较分散，主要分布在流域南部，流域总体生态风险值均较高，低值集聚较少；“低—高”或“高—低”集聚的异常点较少，分布较为规律。说明在研究期间，沱江流域生态风险普遍较高，低生态风险区面积少且分布分散，生态问题形势严峻。

利用软件ArcGIS生成的沱江流域2005—2015年生态风险分布的趋势图(图3)可知，在南北方向上，生态风险在南北方向上呈浅“U”型分布，且北部趋势线上扬，表明生态风险值南北高、中部低，且北部生态风险明显高于南部。在东西方向上，西部趋势线略上扬幅度高于东部，说明西部生态风险略高于中东部。研究区生态风险分布的趋势图和空间自相关集聚的结果(见图2)以及克里金插值的空间分布(见附图5—7)相吻合。

图2 2005-2015年沱江流域生态风险局部自相关特征

注：3个坐标轴X,Y,Z分别代表东西方向，南北方向，生态风险值。图3 2005-2015年沱江流域生态风险空间趋势特征

3.4 生态风险的时空分异

为进一步探究沱江流域在研究期内的空间分异特点，利用软件GS+对生态风险的采样数据进行理论半变异函数的拟合，结果表明2005年球状模型拟合最理想，2010和2015年高斯模型拟合最理想(见表4)。决定系数都大于0.95，说明模型拟合的效果很好。在研究期间，沱江流域生态风险半变异函数的基台值均小于0.01，说明生态风险值与其附近的生态风险值较为集聚，波动较小；2010年的基台值较高，旅游开发建设等人类活动影响较大，2010年之后基台值降低，人们开始追求经济、社会、生态的共同发展，生态风险波动较小。

基于半变异函数的理论模型，运用普通克里金法，对沱江流域在研究期内的生态风险值(ERI)进行插值，得到沱江流域2005—2015年的生态风险空间分布图(附图5—7)，然后统计得出各时期不同等级生态风险的面积(图4)，叠加生成沱江流域研究期间的生态风险转移图(附图8)。

表4 2005，2010和2015年沱江流域生态风险指数半变异函数

图4 2005-2015年沱江流域不同等级生态风险面积变化特征

从不同等级生态风险面积来看，2005年面积占比最大的为较低生态风险，为12 169.48 km2，达到总面积的46.43%；中等生态风险次之，为8 225.37，占比31.38%。相比于2005年，2015年较低生态风险面积由12 169.48 km2下降到了6 040.10 km2，比例下降到了23.05%；中等生态风险面积占比最大，已经超越了总面积的一半，为52.49%；高生态风险区面积增幅较大，由1 896.35 km2增加到了2671.81 km2，增长率达40.89%。总体来说，在研究期间，较低生态风险面积减少，中等生态风险面积增加且增幅较大，较高生态风险面积先增后减，高生态风险面积增加。这是因为在研究期间初始阶段，经济技术水平有限，人们对土地景观的干扰有限，生态环境较少受人类活动的破坏。随着经济的发展和技术的进步，人类活动对土地利用的广度和深度加深，不合理的人类活动对土地利用的结构和生态服务功能产生影响，作为主要生态用地的林地大幅减少，而非生态用地的建设用地大幅增加，在研究期间增幅为96.33%，土地景观趋于破碎化，生态用地的景观连通度降低，生态恢复力降低，生态脆弱度大大提高，景观脆弱性增加，生态风险增加。

就生态风险在行政区划上的分布来说，北部的成都市、德阳市的绝大部分地区在研究期间一直处于高或较高生态风险分布区，因为成都和德阳长久以来就是沱江流域的经济发达地区，工业的发展和郊区农家乐等服务业的兴起，对建设用地需求扩大，建设用地扩张迅速，对单位面积生态服务价值高的林地和耕地产生了侵占，所以生态服务价值总体降低，生态脆弱性增加，生态风险增加。2005—2010年，眉山市和资阳市的大部分地区均为较低生态风险，但是在研究期末期，眉山市和资阳市的生态风险大多转为了中等生态风险；内江市、自贡市、泸州市也由较低生态风险面积占优势转为中等生态风险面积占优势。沱江流域除成都和德阳的生态风险一直处于较高水平，其他地区生态风险水平由低转高，其原因与这些地区经济发展速度的增加是分不开的。

从不同等级生态风险的转移来说，转移面积最多的前三名为：“较低—中等”(较低生态风险转中等生态风险)面积最大，为6 937.81 km2，达到总面积的26.44%；“较高—高”(较高生态风险转高生态风险)的有1 400.08 km2，占总面积的5.34%；“中等—较高”的有636.88 km2，占比2.43%。所有转移中，“较低-低”(较低生态风险转低生态风险)面积最少，只有0.43 km2。由此可知，在研究期间，沱江流域的生态风险转移大多是生态风险低值向高值转移，生态问题严峻。成都市和德阳市在四川省的经济地位突出，二、三产业发达，人口密度大，人为活动对建设用地的需求大，建设用地扩张迅速，景观破碎度高，生态风险等级较高且有向更高级生态风险转移的趋势。

4 结论与讨论

本研究在分析沱江流域土地利用变化的基础上，计算了沱江流域的生态服务价值，在此基础上构建了基于生态服务价值的流域土地景观生态风险评价指标体系，得到沱江流域生态风险指数，并进行半变异函数拟合、克里金插值，定量地分析了研究区基于生态服务价值的生态风险时空分异特征。结果表明：

在研究期内，沱江流域土地利用动态变化显著，主要以建设用地的增加和林地、耕地的减少为主。建设用地由548.94 km2增长到了1 077.75 km2，增幅达96.33%；林地减少了473.60 km2，占林地总面积的11.82%。建设用地面积扩张最多，土地利用转移以耕地转建设用地、林地转耕地为主。这与经济社会的不断发展产生的巨大开发建设需求有关。

沱江流域生态服务价值总量先增后降，由2000年的2.64×1010元增加到了3.93×1010元然后降低到了3.88×1010元，林地、耕地对生态服务价值的贡献最大。在研究期间，林地的贡献率由2000年的62.75%下降到了2015年的37.20%，但一直保持在35%以上，因为林地面积减少，但林地的生态价值系数较高，贡献率仍较高。耕地的贡献率由2.73%增加到了38.43%，这是因为耕地承担的供给功能比较大，虽然耕地面积减少了，但是由于经济的发展，耕作技术的进步，粮食单产提高，耕地对生态服务价值的贡献率一直保持在较高水平。

沱江流域在研究期内生态风险的时空分布有一定规律。从时间上看，随时间的推移生态风险呈逐渐增强的趋势，较低向中等、中等向较高、较高向高生态风险转移；从空间上看，总体上流域北部的成都、德阳的生态风险高，南部的眉山、资阳生态风险较低，自贡、泸州的生态风险逐渐增加。从面积占比来说，生态风险由较低生态风险占优势到以中等生态风险为主。

研究构建了基于生态服务价值的生态风险评价框架，框架由两个方面构成，一是由生态服务价值损失构建的损失度指数，二是由地形、人为胁迫、生态恢复力、景观脆弱性4个维度，坡度、不透水面比例变化、净初级生产力、植被覆盖度、景观连通重要性、景观易损度6个因子构建的概率指标体系得到的概率指数，区别于单纯由景观指数得到的生态风险指数，并且将研究区划分为若干风险小区(格网)，将之与行政区划相叠加，这对各级政府后续的降低生态风险的管控措施的制定具有一定参考价值。