生态系统服务供给是指特定时空范围内生态系统为人类提供的产品与服务, 需求则是指人类消费与享用生态系统生产的产品与服务, 两者将自然生态系统与人类社会系统有机统一起来^[1-2], 成为生态系统保护的重要概念。人类社会的发展在很大程度上通过消费生态系统服务来改善其福祉, 因而生态系统服务的持续供给是区域可持续发展的基础^[3]。当前, 中国正处于快速城镇化阶段, 城市土地的无序扩张还将进一步扩大生态系统服务失衡^[4], 并最终威胁人类福祉和全球可持续发展^[5]。因此, 如何在社会经济发展深度转型期同时兼顾生态系统服务平衡已成为生态经济学、生态学等学科的研究热点^[6-7]。

生态系统服务的定量评估能够反映生态系统功能的强弱, 对区域开展国土空间规划和制定生态补偿政策具有重要指导意义。已有研究对生态系统服务供给进行了大量有益探索, 且计量方法日臻成熟^[8-10]。近年来, 学者们开始关注生态系统服务供给与需求的综合研究^[11-13], 但由于缺乏合理的量化方法, 对于供需平衡关系的研究还尚缺乏。此外, 不少学者认识到城镇化的空间差异是影响生态系统结构、功能和区域差异的重要驱动因素, 并讨论了城镇化对生态系统服务的影响^[14-17]。然而, 对其驱动力的研究忽略了生态系统服务作为一种公共服务具有空间溢出效应^[18-19], 也并未考虑到区域生态系统服务功能的强弱是人为因素和自然因素共同作用的结果^[20-22], 这限制了生态环境保护相关政策的制定和实施。因此, 本研究将在深入探讨两者相互作用关系的基础上, 将自然因素纳入两者关系的分析之中, 使研究结果更具现实意义。

鉴于当前研究的经验与不足, 本文首先采用Burkhard等提出的生态系统服务供需平衡评估矩阵法、城镇化相关指标, 分别核算长江中游地区生态系统服务平衡指数和综合城镇化水平；其次在生态系统服务平衡与城镇化时空特征分析的基础上, 通过双变量空间自相关模型分析两者的空间关系特征；最后利用空间计量模型探索两者的空间依赖、空间溢出效应, 以期为快速城镇化地区生态保护与社会经济协调发展提供理论支持。

1 研究区概况与数据来源 1.1 研究区概况

长江中游地区由湖北、湖南和江西三省组成(108°21′—118°28′E, 20°09′—33°20′N), 总面积为56.47万km², 下辖325个县级行政单元(图 1)。长江中游地区处于中国第二、三级阶梯过渡带, 地形地貌类型复杂多样, 以丘陵山地为主, 属北亚热带季风型气候, 夏热冬温, 四季分明, 雨热同期。在城镇化和工业化快速发展的背景下, 长江中游地区依托“中部崛起”、长江经济带等战略, 逐渐成为我国长江上中下游协调发展的重要链接区和东中西发展的传导区。与此同时, 长江中游地区的土地利用也经历了深刻转型, 1995—2015年间建设用地增加了5852.25 km², 林地减少了1999.53 km², 草地减少了1124.89 km², 耕地减少了4966.64 km²。土地利用变化对生态环境造成了严重干扰, 探索长江中游地区生态系统服务平衡与城镇化的空间关系, 对区域生态环境保护和社会经济可持续发展具有重要意义。

1.2 数据来源

土地利用数据来源于中国科学院资源环境科学与数据中心(http://www.resdc.cn), 空间分辨率为100m。在参照已有土地利用分类体系及区域实际土地利用特点的基础上, 将地类划分为耕地、林地、草地、水域、湿地、建设用地、未利用土地等7个一级类及18个二级类。人口、经济等数据来自各省市统计年鉴。

2 研究方法 2.1 生态系统服务供需平衡测算 2.1.1 土地利用类型划分

由于国内外土地利用分类体系存在差异, 本研究根据国内已有的土地利用分类体系, 对Burkhard等采用的环境信息协调(Co—Ordinated Information on the Environment, CORINE)数据库中的土地利用分类体系进行整合, 以确定本研究的土地分类体系(表 1)。

2.1.2 构建生态系统服务供需平衡评估矩阵

Burkhard等以欧洲生态环境为背景构建了基于专家知识的生态系统服务供需评估矩阵, 该评估矩阵将不同土地利用类型的生态系统服务供给和需求按0到5的区间范围从低到高进行打分, 0代表土地利用类型无生态系统服务供给或需求, 5代表土地利用类型有最高水平的生态系统服务供给或需求。本研究基于Burkhard团队及国内已有研究成果^[23-25], 选取3大类23小类生态系统服务供给/需求值, 利用整合的土地利用类型构建长江中游地区生态系统服务供需评估矩阵(表 2)。

2.1.3 生态系统服务供需平衡指数测算

根据构建的生态系统服务供需评估矩阵, 测算生态系统服务供给指数、需求指数, 进而利用二者的差值表征生态系统服务平衡指数, 各指数计算公式为：

(1)

(2)

(3)

式中, ESSI_t代表生态系统服务供给指数；ESDI_t代表生态系统服务需求指数；ESBI_t代表生态系统服务平衡指数；SM_{ij, t}, DM_{ij, t}, 分别表示在时间t第i种土地利用类型第j种生态系统服务类别的供给、需求值。LUA_{i, t}为时间t第i种土地利用类型的面积。n表示土地利用类型(n=1, 2…, 7), m代表生态系统服务种类。

2.2 城镇化水平测度

城镇化是一个复合系统, 主要包括人口城镇化、经济城镇化、土地城镇化、社会城镇化等四个子系统^[26]。人口增长是人口城镇化的核心要素, 经济增长是经济城镇化的驱动力, 城市土地扩张是空间城镇化的保证, 生活水平的提高是城镇化的最终目标^[16]。由于社会城镇化相对抽象并且难以量化和空间化, 而人口的增长水平、经济发展水平及城市土地比例在一定程度上可以反映社会城镇化, 因此选择人口密度、经济密度、城镇土地比例分别代表人口城镇化、经济城镇化和土地城镇化^[14]。此外, 由于以上三个指标的空间分布高度相似, 因此首先根据公式(4)对三个指标进行标准化, 然后根据公式(5)将这三个指标整合为一个综合指标—综合城镇化水平。

(4)

(5)

式中, PD、GDPD、ULP分别代表人口密度、经济密度、城镇土地比例；x′_ij是第i个研究单元第j个指标的标准值(i=1, 2, 3, …n；j=PD, GDPD, ULP), x_iPD、x_iGDPD、x_iULP分别为第i个研究单元标准化后的指标值, CUL_i是第i个研究单元的综合城镇化水平。

2.3 空间自相关分析

利用双变量Moran′s I检验生态系统服务平衡与城镇化的空间集聚和分散模式。全局双变量Moran′s I用来检验研究区域内生态系统服务平衡与城镇化之间的整体关联程度及其显著性；局部双变量Moran′s I用来识别不同空间位置上可能存在的空间关联模式, 从而把握局部空间要素的集聚和分异特征^[27-28]。二者的计算公式如下：

(6)

(7)

式中, I_G代表全局双变量Moran′s I, I_L代表局部双变量Moran′s I, 取值范围为[-1, 1], 正值表示空间正相关, 负值表示空间负相关, 趋近于0则为空间随机分布；n为空间样本数量；W_ij为基于空间邻接关系建立的空间权重矩阵(ij=1, 2, 3, …n)；z_i^e为第i个研究单元的生态系统服务指数, z_j^u为第j个研究单元的镇化水平。此外, 在局部Moran′s I的LISA聚类图中, 高-高型聚类模式表示高生态系统服务平衡和高城镇化水平值区, 低-低型与之相反, 均为空间正相关; 低-高型聚类模式表示低生态系统服务平衡和高城镇化水平值区, 高-低型则相反, 均为空间负相关。

2.4 空间回归分析

空间溢出效应是指由临近位置引起的空间外部性^[29]。本研究利用空间计量模型探究生态系统服务平衡与城镇化的空间溢出关系。首先利用最小二乘法对二者的空间依赖性进行检验；其次采用包括空间滞后模型和空间误差模型的空间计量模型探究二者的空间关系^[30-31]。两种模型的计算公式如下所示：

(8)

(9)

式中, Y是因变量的向量, X是自变量的N×M矩阵；W是空间权重矩阵, β是回归系数矩阵, μ为空间自回归误差项, ε是空间自回归随机误差项, ρ、λ分别为空间滞后项和空间误差项。

3 结果与分析 3.1 生态系统服务平衡时空特征分析

从时间进程来看, 长江中游地区生态系统服务平衡状况整体呈现出持续恶化的态势, 平衡指数由1995年的47.81降至2015年的46.69, 二十年间共降低了2.40%。利用“自然间断点法”将测算结果划分为5个等级, 以分析各等级生态系统服务供给、需求及平衡的空间分布特征(图 2, 图 3, 图 4)。就空间分布而言, 研究区生态系统服务平衡指数表现为从城市群较为集中的平原向丘陵再向海拔较高的山区逐渐递增的趋势。生态系统服务高平衡区和次高平衡区主要分布在湖北西部的大巴山-巫山-武陵山区, 湖南西部的雪峰山区及南部的南岭山区, 江西东部的武夷山以及江西和湖南省之间的罗霄山区；低平衡区和次低平衡区主要集中在江汉平原、鄱阳湖平原和洞庭湖平原, 尤其是武汉城市圈、环长株潭城市群、环鄱阳湖城市群等三大城市群周围；一般平衡区主要位于山区与平原的过渡地带。生态系统服务供给指数的空间分布与生态系统服务平衡指数的空间分布相似, 而生态系统服务需求指数与生态系统服务平衡指数的空间分布表现出相反的特征。

3.2 城镇化水平时空特征分析

1995—2015年, 长江中游地区城镇化水平表现出明显的空间异质性(图 5)。整体而言, 长江中游地区城镇化水平呈持续上升趋势, 平原地区城镇化水平明显高于周边山区。具体而言, 城镇化水平以武汉、长沙和南昌三大城市为核心向城市群外围递减。此外, 主要交通道路沿线的城镇化水平明显较高, 如京九线、京广线、湘黔线、浙赣线和汉丹线沿线各县, 总体呈“两横三纵”的空间分布特征。

3.3 生态系统服务平衡与城镇化的空间相关性分析

对长江中游地区1995、2005、2015年的生态系统服务平衡与城镇化进行空间自相关分析, 全局Morran′s I指数分别为-0.3221(P=0.001)、-0.3471(P=0.001)、-0.4049(P=0.001), 表现出显著空间负相关性。在Z检验的基础上(P=0.05)绘制双变量局部空间自相关LISA聚集图(图 6), 用于表征区域城镇化水平与生态系统服务平衡之间的空间相关性。1995—2015年, 长江中游地区的主要空间聚类模式为高-低型(高生态系统服务平衡和低城镇化水平), 低-高型(低生态系统服务平衡和高城镇化水平)和低-低型(低生态系统服务平衡和低城镇化水平)。从空间分布来看, 高-低型主要分布在湖北省的大巴山-巫山-武陵山区, 湖南西部的雪峰山区及南部的南岭山区, 江西东部的武夷山以及江西和湖南省之间的罗霄山区, 这些地区森林资源丰富, 生态优势明显。低-高型主要分布在三大平原、城市群及主要交通道路沿线, 这些地区人口密度大, 城镇体系发达, 对生态资源的占用与破坏现象更为严重。低-低型主要分布在周边山区。此外, 该区还有少量高-高型(高生态系统服务平衡和高城镇化水平), 主要分布在长沙市和南昌市周围。

3.4 生态系统服务平衡与城镇化的空间依赖性分析 3.4.1 变量选取

本研究主要探讨城镇化对生态系统服务平衡的影响, 因此空间计量模型中的因变量为1995、2005、2015年生态系统服务平衡指数, 自变量为城镇化水平。此外还选取了一组控制变量以更好地解释城镇化对生态系统服务平衡的影响(表 3)。相关研究证实地形条件在生态系统服务平衡空间格局的形成中起着至关重要的作用^[32]。据此，本研究选取地形位指数代表长江中游地区的地形状况。气候也是影响生态系统平衡的重要空间决定因素, 由于温度和海拔高度之间存在很强的相关性, 因此选择降水来代表气候因素。已有研究表明城镇布局往往倾向于靠近水资源丰富的地带^[33], 因而选择到主要河流水系的距离来表征水文条件。对选取的变量进行方差膨胀因子测试, 结果均小于10, 不存在多重共线性问题。

3.4.2 空间依赖关系分析

表 4显示了1995—2015年长江中游地区生态系统服务平衡与城镇化的空间依赖性统计检验。表中Moran′s I指数在P < 0.05时均显著, 表明空间依赖性存在于所有回归中, 仅采用OLS模型分析会导致结果出现偏差。当在OLS模型中添加控制变量后发现该模型具有更高的R²和对数似然值(Log likelihood)以及更低的AIC和SC值, 模型拟合度的提高表明在研究城镇化对生态系统服务平衡的影响时控制其他影响因素至关重要。此外, 拉格朗日乘子(LM)检测显示LM(滞后)和LM(误差)均较为显著, 表明空间计量模型比OLS模型更能科学的解释二者之间的关系。

空间回归系数表示各选取变量对生态系统服务的影响。结果(表 5)显示1995—2015年城镇化水平与生态系统服务平衡指数为显著负相关关系, 城镇化水平每增加10%, 生态系统服务平衡指数分别减少2.918%, 5.113%, 4.813%。城镇土地扩张过程中不可避免的侵占原有生态土地, 从而导致域生态系统服务失衡。地形位指数与生态系统服务平衡指数为显著正相关关系, 地形位指数每增加10%, 生态系统服务平衡指数相应增加0.314%, 0.501%, 0.398%。随着地形梯度的上升, 人口数量及人类活动逐渐减少, 生态价值较高的林地、草地的比重不断上升, 而建设用地、耕地等生态价值较低的土地类型比重逐渐下降, 致使生态系统服务功能随地形梯度的上升而不断增强。年均降雨量与生态系统服务平衡指数呈显著正相关关系, 年均降雨量每增加10%, 生态系统服务平衡指数相应增加0.236%, 0.567%, 0.225%；降水条件在城镇化对生态系统服务的影响中主要发挥着正向调节作用, 降水丰沛地区的植被拥有较强的自我修复能力，从而有利于植被覆盖的恢复。与其他因素相比, 距河流距离与生态系统服务平衡指数的关联不太显著, 正相关关系表明距离河流越远, 生态系统服务平衡指数越高。这是因为靠近水源是人类定居的首要条件之一, 而距水源较远的土地则受人类活动干扰较少。

4 讨论

(1) 生态系统服务平衡状况整体呈现出持续恶化的态势。

研究期间, 武汉城市圈、长株潭城市群以及鄱阳湖城市群等城市群较为集中的地区对生态系统服务的消耗不断增加, 同时这些地区的生态系统服务供给能力有限, 进而加剧了整个区域的生态系统服务失衡。已有研究表明城镇化会侵占大量生态用地, 从而对生态系统服务产生负面影响^[34-35]。本研究结果显示长江中游地区城镇化水平不断提升的同时生态系统服务平衡状况不断恶化, 再次印证了该结论。

今后, 应继续大力推进长江流域生态文明建设, 严格实施退耕还林、还湖以及湿地保护等生态政策, 努力实现湿地、水域不萎缩, 功能不衰减, 生态不恶化。

(2) 不同区域生态系统服务对城镇化进程的响应具有空间异质性。

双变量局部自相关分析发现, 除高-低型、低-高型两种主要聚类模式外, 研究区还有高-高型、低-低型聚类模式。高-低型、低-高型聚类模式表示生态系统服务与城镇化之间的负相关关系, 而低-低型、高-高型聚类模式则表明区域内的地形条件、植被覆盖率、降水量等其他非城镇因素对区域生态系统服务调节也起到重要作用^[36]。例如南昌县境内赣江流域沿线高楼林立、人口稠密, 但由于水系和湿地覆盖率较高, 区域在生物多样性、户外文化娱乐等方面具有较高的服务价值^[37], 致使该县为高-高型城镇。又例如修水县的城镇化水平较低, 同时生态系统服务平衡也较低, 这与该县较高的地形坡度和较强的降水有着紧密联系。已有研究指出区域自然资源本底以及社会经济条件的差异会导致生态系统服务对城镇化的响应有所不同^[38], 与本文结论一致。

基于此, 本文建议各区域有必要因地制宜探索城镇化建设道路。高-低型城镇在继续充实基础设施的同时应大力发展生态旅游、特色农产品加工等绿色产业, 努力将自然生态优势转化为经济优势。低-高型城镇应继续保持经济发展速度, 同时控制人口规模、减缓城市扩张速度, 以供给侧结构性改革为主线助推经济社会持续健康发展。基于“生态保护优先”原则, 应考虑将高-高型城镇划为“生态经济区”，在这些区域应鼓励新兴产业、高新技术产业发展, 强制淘汰污染型产业。低-低型城镇在继续推进城镇化的同时, 还应在生态退化较为集中的区域布置生态修复工程, 如植树造林、边坡绿化等，努力实现经济与生态环境的协调发展。

(3) 生态系统服务平衡与城镇化水平存在空间溢出效应。

本研究将城镇化水平及各自然因素与生态系统服务平衡进行空间回归分析得出, 生态系统服务平衡不仅受到研究单元自身因素的影响, 而且还会受到临近单元相同因素的干扰。总体而言, 城镇化对生态系统服务的影响表现为空间负外部性^[39]。这种负外部性可能是城镇化地区与临近区域持续进行物质和能量交换的结果。当某一城镇化地区的生态环境要素(如地形、气温、降水、二氧化碳等)发生变化时, 将通过大气循环、水循环以及生物迁徙等方式影响临近区域的生态系统服务平衡状况^[40]。

基于该结论, 本研究建议各研究单元在制定区域发展规划时充分考虑城镇化的负外部效应, 避免在与自然生态系统直接相邻的地区进行生态破坏活动。同时, 尝试加强县际合作, 实现生态与经济资源的共享与互补。

(4) 后续待开展的工作。

本研究基于县级尺度从宏观上探讨了长江中游地区生态系统服务平衡与城镇化的空间关系, 研究结果对于区域制定总体发展规划和生态保护相关政策具有一定借鉴意义。然而, 各县级行政单元内生态系统服务平衡与城镇化的相互关系还不明晰, 不利于县级相关规划的实施与落地, 因此下一步工作拟选取典型县市从格网尺度深入分析二者之间的内在关系。

5 结论

本研究以城镇化对生态系统服务平衡的影响为出发点, 采用生态系统服务供需平衡评估矩阵法、双变量空间自相关分析模型以及空间计量模型等技术方法, 探索了生态系统服务平衡与城镇化水平的相互作用关系, 研究结论如下：

(1) 研究期间长江中游地区生态系统服务平衡状况不断恶化, 周边山区是生态系统服务平衡高值区, 腹部三大平原为生态系统服务平衡低值区。

(2) 研究期间长江中游地区城镇化水平呈上升态势, 地处江汉平原、洞庭湖平原、鄱阳湖平原的武汉市、长沙市、南昌市等三大核心城市及主要交通道路沿线是城镇化水平的高值区, 周边的雪峰山、南岭和武夷山等山区为城镇化水平的低值区。

(3) 研究区生态系统服务平衡与城镇化水平之间存在空间依赖关系, 且不同区域的生态系统服务对城镇化的响应有所差异, 生态系统服务平衡较低的区域更容易受到城镇化的侵扰。

(4) 添加控制变量的空间计量模型更好的解释了城镇化对生态系统服务平衡的影响, 结果显示城镇化对态系统服务平衡具有空间负外部效应, 这种负外部效应在城镇化发展到一定阶段时逐渐减弱。