随着全球化的不断深入以及全球范围竞争的不断加剧, 城市群日益成为全球经济竞争的新空间单元^[1-3]。全球范围内, 城市群地区是区域经济、城镇化和工业化快速发展的地区, 同样也是社会问题、资源枯竭、环境恶化和土地利用问题等一系列矛盾最为突出的地方^[4-6]。城市群的发展不仅对区域经济发展至关重要, 对社会进步以及提高国家国际影响力也十分重要^[2]。另外当前中国政府将城市群建设作为政策工具在全球和地区竞争中占据有利地位, 并且在未来十年中制定了数十个国家级、地区和地方级城市群计划, 但是对于城市群的研究中空间视角往往被忽略^[3], 城市群地区快速的城镇化发展带来的土地利用/土地覆被变化的生态环境效应研究仍然不全面, 以单一城市群生态系统为研究对象不能从根本上解决区域生态问题, 无法制定城市圈集聚区地区合理的生态环境和土地利用政策^[6]。因此，很有必要研究城市群集聚地区城镇化发展对生态系统服务的影响，这对于区域可持续发展具有重要意义。

城市群的发展必然会导致土地利用方式、格局和程度的变化, 对区域生态系统构成严重威胁, 因此对城市群地区生态系统服务的影响进行建模分析就显得十分必要^[7]。快速城镇化进程中的土地利用变化及其生态环境效应已成为生态系统服务共同国际分类(Common International Classification of Ecosystem Services, CICES)、千年生态系统评估(Millennium Ecosystem Assessment, MEA)、生态系统和生物多样性经济学(Economics of Ecosystems and Biodiversity, TEEB)、生物多样性和生态系统服务政府间科学政策平台(Intergovernmental Science-Policy Platform on Biodiversity and Ecosystem Services, IPBES)和国际土地变化科学(Land Change Science, LCS)及国际科学理事会(International Council for Science, ICSU)最新大型科学计划“未来地球”(Future Earth)等组织和倡议的重点关注内容和研究主题^[7-8]。在此背景下, 探寻城市群集聚区土地利用/土地覆被变化的生态环境效应的空间溢出特征对于跨区域土地利用规划和生态保护协同治理具有重要意义。

生态系统服务是指人类直接或间接从生态系统得到的利益, 主要包括供给服务(如提供食物和原材料)、调节服务(如水调节和气候调节)、文化服务(如娱乐和休闲文化)以及支持服务(如维持生物多样性和保持土壤)^[8-9]。生态系统服务关乎人类福祉, 是人类生存和经济社会发展的基础^[9-10]。生态系统服务长期以服务流的形式出现, 能够带来这些服务流的是自然资本, 这种服务流可以增进人类的福祉^[9]。城市群地区快速的城镇化和工业化导致了土地利用的快速转型^[11-13], 严重的干扰了生态系统的结构和功能^[14-15]。据Millennium Ecosystem Assessment(2005)研究结果, 在过去的几十年中, 世界上大约60%的陆地生态系统服务功能在退化, 未来还将要持续恶化^[10]。全球范围对于快速城镇化下生态系统服务进行了广泛的评估^{[14, 16-19]}, 目前生态系统服务主要通过价值量^{[9, 20]}、能值^[21]和物质量^[22-25]等三方面进行评估。三种评估方法各有优缺, 但是都能为决策者和管理者进行决策和规划提供关键和有效信息。其中价值量评估方法在评估过程中往往存在主观性强、结果具有不确定性的问题, 但是价值量评估法便于各项生态系统服务之间比较, 评估结果可纳入国民经济核算体系, 进而为生态环境保护、生态功能区划、环境经济核算和生态补偿决策提供重要依据和基础^[26]；另外价值量评估可以有效体现生态系统服务价值的稀缺性, 为国际合作和政府决策提供科学参考^{[9, 27]}。

城市群是城市发展到成熟阶段的最高空间组织形式, 城市群逐渐成为新型城镇化的主体形态和现代化建设的重要载体, 探索城市群地区城镇化与生态系统服务的空间关系对于城市群地区生态保护和土地可持续利用提供科学的借鉴。以往研究对城市群地区开展了大量的生态系统服务评估、生态系统服务权衡、生态系统服务影响因素、生态系统服务供需分析^[28-30], 但是仍存在很多不足, 包括尺度效应、空间相互作用、空间依赖效应和空间溢出效应, 这限制了生态保护政策和土地利用政策的有效制定和实施。现有研究大都基于行政单元的统计数据, 往往掩盖生态系统内部空间异质性特征以及人口分布和经济发展的空间差异, 本研究使用空间化数据可以开展格网尺度上的定量研究, 更能凸显其现实意义。

基于以往研究本文从三方面进行了创新：(1)研究尺度方面, 本文选择长江中游城市群地区作为研究区, 长江中游城市群是一个城市簇, 对于城市群集聚区研究可以更好地为跨区域的协同治理提供科学指导。另外本文采用ArcGIS 10.3中Fishnet工具生成5 km×5 km和10 km×10 km的格网, 作为测度生态系统服务和城镇化的水平基本单元, 探测生态系统服务和城镇化之间的空间相互作用。(2)本研究从格网尺度评估了长江中游城市群地区的城镇化发展水平, 打破以行政单元为评价单元的传统。(3)借助双变量空间自相关模型(Bivariate Moran′s I)、普通最小二乘法(Ordinary Least Square, OLS)、空间滞后模型(Spatial Lag Model, SLM)、空间误差模型(Spatial Error Model, SEM)和加入空间滞后项的空间误差模型(Spatial Error Model with Lag Dependence, SEMLD)模型对长江中游城市群地区城镇化水平和生态系统服务价值之间的空间相互作用、空间依赖效应和空间溢出效应进行了测度。

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

长江中游城市群是中国的国家城市群之一, 位于其他4个国家城市群的中心(长三角城市群、珠三角城市群、成渝城市群和京津冀城市群)。长江中游城市群包括武汉城市群、宜荆荆城市群、长株潭城市群和环鄱阳湖城市群(26°07′05″—30°23′N, 110°15′—118°28′58″E)(图 1)。长江中游城市群属于亚热带季风性气候, 雨量充沛、日照充足、四季分明。地处长江中游, 江汉平原、洞庭湖平原和鄱阳湖平原分布其间, 土壤肥沃, 是中国重要的粮食产地。随着中部崛起战略和长江中游经济带战略的实施, 长江中游城市群已经成为中国经济发展新增长极。另外, 长江中游城市群作为中西部新型城镇化先行区、中国低碳经济发展先行区、内陆开放合作示范区和"两型"社会建设引领区, 在城镇化建设以及生态环境保护方面处于领先地位。快速的城镇化发展背景下, 长江中游城市群土地利用发生了深刻转型, 带来了显著的生态环境变化。鉴于此, 选择长江中游城市群作为研究区, 探索长江中游地区土地利用变化生态环境效应与城镇化水平之间的空间关系, 对于城市群集聚地区土地可持续利用和生态环境保护具有重要的借鉴和指导意义。

1.2 数据来源

研究涉及1995年、2005年和2015年土地利用/土地覆被变化数据来源于源于中国科学院资源环境科学数据中心(http://www.resdc.cn), 空间分辨率为30 m×30 m^[31-32]。借助ArcGIS 10.3工具, 将二级地类重分类为耕地、建设用地、林地、草地、水体、湿地和未利用地等7种土地利用类型。1995年、2005年和2015年1 km分辨率人口和GDP栅格数据来源于中国科学院资源环境科学数据中心(http://www.resdc.cn)。粮食单产数据来源于1996—2016年份《湖北省统计年鉴》、《湖南省统计年鉴》和《江西省统计年鉴》, 粮食价格数据来源于2016年《中国农产品价格调查年鉴》等统计资料。基础地理数据来源于国家基础地理信息系统全国1:400万数据库(http://nfgis.nsdi.gov.cn)。

2 研究方法 2.1 生态系统服务价值测度

Costanza等提出的生态系统服务分类方法以及当量表不适于直接用于中国生态系统服务评价, 可能会低估或忽视某些类型生态系统服务^{[9, 20]}。谢高地等基于700多位生态学者的意见, 对中国生态系统服务分类和当量表等进行了修正^[20]。1个标准单位生态系统服务价值当量因子定义为1 hm²农田每年平均粮食产量的经济价值的1/7。其他生态系统类型生态系统服务价值当量确定依据该类生态系统相对于农田生态系统的重要性。由于土地利用类型不能和生态系统类型一一对应, 本研究选择最接近的土地利用类型进行当量赋值。本研究参考陈万旭等^[29]对长江中游地区的研究结果, 计算出生态系统服务价值当量因子为2148.35元hm^-2 a^-1。公式如下：

(1)

(2)

式中, ESV代表生态系统服务价值；AESV代表地均生态系统服务价值；LUC_i代表i类土地利用类型面积；VC_ij是第i类土地利用类型j类的生态系统服务当量值。m是生态系统服务类型数量(m=9), n是研究区土地利用类型的数量。

2.2 城镇化水平测度

城镇化是一个复杂的系统, 主要包括人口城镇化子系统、经济城镇化子系统、社会城镇化子系统和土地城镇化子系统等四方面, 任何一方面的偏颇都将降低城镇化综合效益^[33-35]。其中人口城镇化是城镇化的核心, 经济城镇化是动力, 土地城镇化是保障^[33]。人口城镇化和经济城镇化为城镇化和生活方式的转变奠定基础, 社会城镇化和土地利城镇化分别是社会和空间的表现形式^{[15, 36]}。由于社会城镇化社会经济数据获取比较困难, 本文不与讨论。本文选择人口密度、经济密度和建设用地面积比重格网数据来量化人口城镇化、经济城镇化和空间城镇化。由于这3个指标的高度相关性, 本文采用极差标准化方法对这3个指标进行标准化, 然后进行等权重加权计算得到城镇化发展水平。

2.3 空间自相关分析

探索性空间数据分析方法的核心内容是探索某一属性值是否在空间上表现出分布的聚集性和空间异常等特征, 一般通过全局Moran′s I指数和局域Moran′s I指数(LISA)对空间关联模式(趋同或异质)进行度量和检验, 进而揭示对象的空间分布特征^[37]。全局空间自相关用来验证区域某种属性值在整个区域空间分布的集聚态势, 常用Global Moran′s I指数来测度^[38]。公式如下：

(3)

式中, x_i, x_j为观测值；为x_i的平均值；W_ij为空间单元i与j的空间权重邻接矩阵(ij=1, 2, 3, …n)。全局Moran′s I取值一般在-1到1之间, 若大于0, 说明空间存在正自相关；若小于0, 说明空间实体呈离散分布；若其数值趋近0, 则说明该空间服从随机分布, 常用P值进行显著性检验。

本文基于Anselin和Rey^[37]提出了基于双变量的区域空间自相关探索城镇化和生态系统服务价值之间的空间关系。双变量空间自相关分析所产生的Moran′s I值是用所有相邻位置的加权平均值评估一个位置变量值与其他变量的相关程度^[38-39]。公式如下：

(4)

式中, w_ij为空间单元i、j之间的空间连接矩阵、；x_kⁱ是空间单元i属性k的值；x_lⁱ是空间单元j属性l的值；、是属性k、l的平均值；λ_k、λ_l是属性k、l的方差。根据Local Moran′s I指数, 计算结果中H-H类型表示高地均生态系统服务价值和高城镇化水平类型区, L-L类型表示低地均生态系统服务价值和低城镇化水平类型区, H-L类型表示高地均生态系统服务价值和低城镇化水平类型区, L-H类型表示低地均生态系统服务价值和高城镇化水平类型区。

2.4 空间回归分析

以往研究中很少考虑到生态系统服务作为一种公共服务具有很强的空间溢出效应^{[8, 40]}, 空间生态系统服务空间无关联性及均质性的假设会导致研究结果的不完整和不科学^[41]。本文引入空间经济计量模型来解决生态系统服务价值和城镇化之间的复杂空间相互作用与空间依赖性问题^[42-43]。首先采用普通最小二乘法从整体上探测生态系统服务价值和城镇化之间的关系, 而不考虑邻近区域的空间影响。然后再采用适用于截面数据的空间常系数回归模型, 包括空间滞后模型、空间误差模型和加入空间滞后项的空间误差模型三种空间回归模型探索地均生态系统服务价值和城镇化水平之间的空间关系^{[30, 44-46]}。OLS模型和空间回归模型均在GeoDa 095i中实现。

3 结果与分析 3.1 长江中游城市群生态系统服务价值时空特征

研究表明1995—2015年间长江中游城市群生态系统提供的生态系统服务价值总体呈现降低趋势, 由1995年的13505.82亿元降低到2015年的13513.86亿元(表 1), 其中1995—2005年增加了0.40%, 2005—2015年降低了0.33%。研究期间耕地提供的生态系统服务价值不断降低, 研究期间降低了93.14亿元；林地和草地提供的生态系统服务价值同样呈现持续降低趋势, 而且降低幅度持续增加；1995—2000和2000—2005年间林地提供的生态系统服务分别降低了0.20%和0.82%, 草地提供的生态系统服务分别降低了2.53%和4.61%。1995—2015年间长江中游城市群生态系统服务价值增加的主要原因是水体和湿地提供的生态系统服务价值的增加。未利用地对长江中游城市群生态系统服务价值贡献最低。长江中游城市群生态系统功能结构中水文调节功能占比最大(18%左右), 其次是维持生物多样性功能(14%左右), 食物生产功能占比最小(3%以下)(表 2)。长江中游城市群地区林地比重较大, 为水文调节提供了可能, 另外长江中游城市群地区地形复杂, 水体和湿地资源丰富, 为生物资源多样性提供了可能。快速的城镇化占用了大量的耕地是食物生产功能降低的主要原因。研究期间其中食物生产功能、原材料生产功能、气体调节功能和保持土壤功能持续降低, 1995—2015年间分别降低了2.90%、1.24%、1.20%和1.68%。提供美学景观、废物处理和水文调节功能在研究期间有显著的增加。

从地均生态系统服务价值空间分布来看, 地均生态系统服务价值高值区主要集中在长江沿线地区, 在鄱阳湖和洞庭湖地区尤为显著, 另外在中部的幕阜山、罗霄山区以及周边山区也是生态系统服务价值高值区；低值区主要分布在江汉平原、洞庭湖平原和鄱阳湖平原地区以及以武汉市、长沙市、南昌市以及宜昌市等地级市以及周边格网。从生态系统服务价值变化率来看, 主要的交通道路沿线地区是生态系统服务价值恶化的主要区域, 另外主要的地级市以及周边格网同样是生态系统服务价值降低较为显著的地区(图 2、图 3)。主要是由于这些区域是城镇化水平相对较高, 同时对生态系统的干扰也较为强烈。就不同尺度的生态系统服务价值变化来说, 1995—2005和2005—2015年间5 km尺度生态系统服务价值降低格网数量比重分别为45.64%和55.62%, 10 km尺度生态系统服务价值降低格网数量比重分别为44.27%和59.08%。可以看出生态系统服务价值具有显著的尺度效应, 另外长江中游地区生态系统服务恶化面积在逐渐增加。

3.2 长江中游城市群城镇化水平时空特征分析

总体来说, 长江中游城市群地区山区城镇化水平低于平原地区, 具体来说城镇化水平以武汉市、长沙市、南昌市和宜昌市等地级市为核心向周边逐渐减弱。就城镇化水平的变化而言, 1995—2005年和2005—2015年间长江中游城市群地区城镇化水平显著增加区域呈现面状和线状分布, 大城市周边呈现面状分布, 主要的道路交通沿线地区呈现线状分布(图 4、图 5)。对比不同尺度的城镇化水平变化来说, 1995—2005年和2005—2015年间5 km尺度城镇化水平增加格网面积比重分别为81.76%和97.34%, 10 km尺度城镇化水平增加格网面积比重分别为81.69%和98.69%；同样可以看出城镇化水平具有显著的尺度效应, 长江中游地区城镇化水平增加格网面积总体上在逐渐增加。

3.3 长江中游城市群生态系统服务价值与城镇化之间的空间相关性特征

基于GeoDa 095i软件, 测度了1995年、2005年和2015年间不同格网尺度生态系统服务价值与城镇化水平之间的双变量全局空间自相关Moran′s I指数。1995年、2005年和2015年5 km格网尺度生态系统服务价值与城镇化水平之间的双变量全局空间自相关Moran′s I指数分别为-0.274、-0.264和-0.275, 10 km格网尺度生态系统服务价值与城镇化水平之间的双变量全局空间自相关Moran′s I指数分别为-0.293、-0.271和-0.263, 而且通过了0.001水平的显著性检验。研究期间城镇化水平和生态系统服务价值之间具有显著的空间负相关, 表明城镇化的发展会导致生态系统服务的恶化。对比不同格网尺度, 1995年和2005年5 km尺度生态系统服务价值变化和城镇化水平变化之间全局Moran′s I指数显著低于10 km尺度, 而2015年则高于10 km尺度。双变量空间自相关检验表明生态系统服务价值和城镇化之间具有显著的空间依赖效应, 研究期间产生了显著负向的外部性。主要是由于随着城镇化的推进, 期建设用地迅速扩张, 对生态系统服务造成强烈的干扰。空间溢出效应的存在主要是由于物质流、能量流和信息流在相邻区域之间的流动而形成。另外自然过程(例如, 大气循环、气温、降水)和人类活动(例如, 城镇化、经济发展、人口迁移)同样会导致生态系统服务的空间溢出。因此生态系统的保护以及治理不能局限于一个单元, 跨区域以及不同等级政府部门的协同治理是生态系统保护的有效手段。

图 6和图 7是1995年、2005年和2015年间长江中游城市群5 km和10 km尺度生态系统服务价值和城镇化水平之间的双变量局部空间自相关LISA聚集图。对比不同时间点5 km和10 km尺度生态系统服务价值和城镇化水平之间的双变量局部空间自相关LISA聚集图(图 6、图 7)可以发现一些类似的空间特征。研究期间, H-L类型(高地均生态系统服务价值和低城镇化水平类型区)和L-H类型(低地均生态系统服务价值和高城镇化水平类型区)是主要的关系类型。其中H-L类型主要分布在周边山区以及中部地区的幕阜山、罗霄山区, 而L-H类型主要分布在平原地区的大城市周边地区(例如, 武汉市、长沙市和南昌市)。H-H类型(高地均生态系统服务价值和高城镇化水平类型区)和L-L类型(低地均生态系统服务价值和低城镇化水平类型区)在长江中游城市群地区也有出现, 但是数量不多。城镇化过程中人口和经济的集聚必然会带来建设用地的扩张, 城市中建设用地的扩张会带来生态系统服务的显著降低, 但是城市中水体、湿地以及绿地可以在很大程度上抵御外界对生态系统的干扰, 城镇化和生态系统保护的协调发展, 可以显著降低城镇化对生态系统服务的干扰。

3.4 长江中游城市群城镇化水平对生态系统服务价值的影响分析

基于生态系统服务价值与城镇化水平之间的双变量全局空间自相关性分析可以看出生态系统服务和城镇化水平之间存在显著的空间依赖性。为了进一步研究二者的关系, 本研究采用了一系列的回归模型来研究城镇化水平对生态系统服务的影响。首先在不考虑空间相互作用的情况下, 进行OLS分析。OLS回归分析结果表明生态系统服务价值变化和城镇化水平变化之间存在显著的空间依赖性, 因此采用OLS回归模型会导致得出的结果不能完整、科学的解释两者之间的关系。另外加入城镇化水平平方项的模型表现优于不加入平方项的模型, 因此加入城镇化水平的平方项能更好的揭示二者之间的相互作用关系(表 3)。空间回归结果表明空间回归模型显著地提高了模型解释能力, 本研究应该采用空间回归模型。对比SEM、SLM和SEMLD三种空间回归模型, 可以发现SEMLD模型具有更大的Log likelihood值, 更小的SC和AIC值, 能够更好地解释生态系统服务价值和城镇化水平之间的空间关系(表 4和表 5)。从回归结果来在所有的模型中都可以发现生态系统服务价值和城镇化水平之间存在显著的负相关关系, 都在0.001水平上显著。可以看出城镇化水平增加对生态系统服务价值具有负向的影响, 即城镇化水平的增加会带来生态系统服务价值的降低。用SEMLD模型来解释说明, 1995年、2005年和2015年5 km尺度城镇化水平每增加1%会导致生态系统服务价值降低0.481%、0.351%和0.363%, 10 km尺度城镇化水平每增加1%会导致生态系统服务价值降低0.547%、0.475%和0.445%。另外可以发现, 本研究中所有的城镇化平方项都在0.001水平上显著, 表明城镇化和生态系统服务之间存在U型曲线曲线, 也就是说在城镇化的初期, 城镇化的增加会导致生态系统服务的降低, 而在城镇化水平达到一定的阈值后, 城镇化的增加会促进生态系统服务的增加。对比不同格网尺度的回归系数可以看出, 5 km格网尺度的回归系数显著低于10 km格网尺度回归系数。所有的模型中的Spatial lag term都在1%水平上显著, 表明生态系统服务价值变化不仅受到自身格网要素的影响, 而且受到相邻或者更远格网单元的要素的影响。1995年、2005年和2015年5 km格网尺度周边地区生态系统服务价值平均每增加1%会导致自身格网单元生态系统服务价值增加1.034%、1.044%和1.033%, 10 km尺度周边地区生态系统服务价值平均每增加1%会导致自身格网单元生态系统服务价值增加1.030%、1.032%和1.021%。另外, 所有模型中Spatial error term同样在1%水平上显著, 表明生态系统服务价值变化也受到其他非城镇化因素(地形、气候、交通等)的影响。

4 结论与讨论 4.1 结论

基于长江中游城市群1995年、2005年和2015年3期土地利用遥感监测数据, 借助生态系统服务价值评估方法、双变量空间自相关分析以及空间回归分析模型等技术分析手段, 分析并揭示了长江中游城市群地区生态系统服务价值和城镇化之间的空间关系。研究结果如下：

(1) 研究期间长江中游城市群地区生态系统服务价值从由1995年的13505.82亿元降低到2015年的13513.86亿元。长江沿线水域以及中部的幕阜山、罗霄山区以及周边山区是生态系统服务价值高值区, 江汉平原、鄱阳湖平原和洞庭湖平原以及主要的地级市周边地区是生态系统服务价值的低值区。

(2) 研究期间80%以上的格网城镇化水平在增加, 山区城镇化水平远低于平原地区, 城镇化水平以武汉市、长沙市、南昌市和宜昌市等地级市为核心向周边逐渐扩散。

(3) 研究期间地均生态系统服务价值和城镇化水平之间具有显著的空间依赖性, 总体上来说城镇化的发展会导致生态系统服务的降低。但是城镇化对生态系统服务的影响存在阈值, 城镇化达到阈值之前城镇化的发展会导致生态系统服务的降低, 而在城镇化发展水平超过阈值以后, 城镇化的发展会促进生态系统服务的改善。

(4) 空间回归模型对城镇化和生态系统服务之间具有较好的解释能力, 回归结果显示城镇化水平的增加会导致生态系统服务价值的降低, 5 km格网尺度的回归系数显著低于10 km格网尺度回归系数。

4.2 讨论

(1) 本研究从5 km和10 km格网尺度来测度长江中游城市群地区的生态系统服务价值和城镇化发展水平之间的相互作用关系, 打破了以往以行政区为研究单元的传统。生态系统服务价值和城镇化空间分布特征以及二者之间的空间关系具有显著尺度效应。

(2) 城镇化过程伴随着人口和经济的集聚以及城市建设用地的扩张，城市建设用地的扩张必然会导致生态系统服务的降低，但是人口和经济的集聚不一定会导致生态服务的退化。这主要是由于人口和经济的集聚不能直接对生态服务产生干扰，而是通过促进建设用地的扩张影响生态系统服务。城市建设用地扩张会降低生态服务，但是城市内部水域、湿地和绿地会在一定程度缓解生态系统功能的退化。

(3) 生态系统服务价值不仅受到自身格网单元要素的影响, 而且受到邻近格网单元或者是更远格网单元要素的影响。因此在未来土地利用规划和生态系统保护规划过程中, 需要综合考虑自身网单元要素以及邻近格网单元或者是更远格网单元要素的影响, 跨区域的协同治理会提高土地利用规划和生态保护规划的实施效果。