人类活动对陆地表层系统的改造导致生态环境的变化, 继而对生态系统服务产生影响^[1], 而建设用地是人类活动作用于陆地表层系统最直接的体现^[2]。城乡建设用地是城乡发展的空间载体, 快速城镇化和工业化推动了城乡建设用地大规模扩张, 从而引发耕地流失、景观破坏、环境污染等生态环境问题^[3—4]。城乡建设用地关联是指城镇建设用地和农村居民点用地在时空维度上的相互联系, 强调从整体上把握城乡建设用地的格局与特征。国土空间治理体系的变革倒逼粗放低效式的建设用地扩张模式加快转型^[5—6], 如何定量刻画城乡建设用地的关联特征, 并探究其对生态系统服务的影响, 是新时期打造高质量国土空间格局的关键问题。

关于城乡建设用地的相关问题备受国内外学者的关注^{[2, 7]}, 其中国外相关研究集中于建设用地方面^[8—11], 尤其是建设用地扩张与人口、经济、政策等要素的相互关系^{[9, 11]}；国内相关研究集中于城乡建设用地的时空演变^[12]、驱动机制^[13—14]、模拟预测^[15]、转型发展^[16]等方面。而随着建设用地扩张带来的负面效应越来越明显, 人地矛盾愈发突出, 关于城乡建设用地利用效率^[17—18]、与人口经济的协调关系^[19—20]、生态环境效应^[21—22]等问题的研究逐渐丰富。其中, 城乡建设用地扩张对生态环境的影响尤其被重视, 不少学者从建设用地对生态、农业空间挤占的时空变化视角出发, 探析其对植被净初级生产力、气候调节、粮食生产、生物多样性等方面的影响^[23—26], 进而分析其对生态系统服务的影响^[27—28]。研究方法上, 利用地理信息技术, 对城乡建设用地和生态系统服务实现时空分析^[29—31], 使得研究更加直观, 并且能够从非线性关系上探讨两者之间的空间关系。整体来看, 现有研究多从单一视角对城镇建设用地或农村居民点用地进行分析, 缺乏关联视角下城乡建设用地时空特征的探讨；同时, 分析城乡建设用地扩张对生态系统服务非线性影响的研究较为少见。

洞庭湖地区是典型的水陆交错生态脆弱带^[32], 也是我国重要的粮食生产基地。城乡建设用地扩张导致土地利用功能发生转型, 给粮食生产和生态安全带来了巨大挑战。如何准确识别城乡建设用地关联的时空特征及其生态系统服务效应, 促进城镇化发展、粮食保障和生态保护统筹协调, 是其区域城乡空间发展中面临的关键问题。因此, 以洞庭湖地区为研究对象, 基于时空关联视角运用引力模型、双变量空间自相关、面板数据模型等方法, 定量识别城乡建设用地关联的时空特征, 并揭示其对生态系统服务的影响, 以期丰富关联视角下土地利用研究的理论框架, 为洞庭湖地区的国土空间高质量发展提供理论支撑。

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

洞庭湖地区位于湖南省北部, 是长江中下游重要的调蓄性湖泊生态区和农业主产区, 也是长江经济带的重要组成部分。其以洞庭湖为中心, 向周边依次过渡为平原、台地和低山丘陵^[33], 大部分地区海拔在30—50m之间, 属亚热带季风气候, 区内河网密布, 水、气、土、生物资源要素丰富, 且处于长江水系共轭与分异的交错并存之节点, 具有显著的生态脆弱性和敏感性, 经济社会发展与生态环境承载之间存在突出矛盾^[34], 探究其城乡建设用地关联的生态环境效应具有重要的现实意义。本文以环洞庭湖地区岳阳、常德、益阳三市的市辖区(岳阳楼区、君山区、云溪区、武陵区、鼎城区、赫山区、资阳区)和岳阳县、华容县、湘阴县、汨罗市、临湘市、汉寿县、安乡县、澧县、临澧县、津市市、南县、沅江市为研究区域, 共计15个县(区、市), 如图 1。

1.2 数据来源

土地利用数据：不同土地利用类型的食物产量、碳密度和生境质量存在差异, 因此为了更加精确计算研究区范围内的食物产量、碳密度和生境质量, 将土地利用类型分为：农田、灌木林、常绿阔叶林、落叶阔叶林、常绿针叶林、落叶针叶林、混交林、草地、疏林地、沼泽地、建设用地、裸地、水体。空间分辨率为0.3 km。

食物产量(粮食作物、油料作物、糖料作物、肉类、水果、牛奶)数据：来源于《中国统计年鉴(1996—2019)》和《中国农村统计年鉴(1996—2019)》；各类食物热量成分数据源于美国农业数据库和相关研究^[35], 见表 1。

碳密度数据：根据不同研究结论得到不同碳库的碳密度数据^[35]。需要说明的是, 由于不同地类地下死亡有机体碳密度数据不全, 因此本文没有将地下死亡有机体碳库纳入碳储量计算范围。不同土地利用类型碳密度, 见表 2。

生境质量数据：水田、旱地、建设用地、裸地是人类活动频繁或自然环境较为恶劣的地类, 因此在本文中将该四类用地作为威胁源。各地类的生境适宜性、胁迫因子权重以及最大胁迫距离、生境类型对胁迫因子的敏感性等相关参数设置来自相关文献^[35], 见表 3和表 4。

控制变量：归一化差异植被指数从地理空间数据云(http://www.gscloud.cn/)获取, 分辨率为1km；夜间灯光、人口密度、GDP密度和交通网络数据从中国科学院地理科学与资源研究所资源环境数据云平台(http://www.resdc.cn/)获取, 其中夜间灯光、人口密度和GDP密度的分辨率为1 km。

1.3 研究方法 1.3.1 引力模型

本文基于改进的引力模型建立城乡建设用地关联模型^[36—37], 以测算城乡建设用地关联强度, 计算公式如下：

(1)

式中：F_u-r代表城乡建设用地时间关联强度；ULP代表网格的建设用地比例；GDPD代表网格的GDP密度；POPD代表网格的人口密度；Dis_u-r代表农村居民点用地网格到城镇建设用地网格的距离。

1.3.2 双变量空间自相关模型

本文利用双变量空间自相关模型(Moran′s I)描述城镇建设用地与农村居民点用地的空间关联, 包括全局Moran′s I和局部Moran′s I。全局Moran′s I用于检验城镇建设用地与农村居民点用地之间的空间相关性是否存在, 局部Moran′s I用于检验不同的网格单元是否存在空间相关性^[38]。计算公式如下：

(2)

(3)

式中：I_u-r和I′_u-r分别代表城乡建设用地全局双变量Moran′s I和局部双变量Moran′s I；W_ij是一个n×n的空间权重矩阵；z_i^u是第i个网格的城镇建设用地占比值；z_j^r是第j个网格的农村居民点用地占比值。

1.3.3 生态系统服务测度方法

生态系统服务(ESs)是生态环境或生态系统与人类活动进行交互并发挥其内在功能作用的主要途径和方式, 其常常通过生态系统服务价值(ESV)来表征, 也即通过ESV来表示ESs的质和量水平^[39]。城市化和工业化进程改变了生态系统的结构(如林地、耕地等)及过程(如生境质量), 使自然生态系统转变为以人为主导或人与自然相耦合的生态系统, 对洞庭湖地区的可持续发展产生了显著影响。针对这些问题, 本文选择了三种典型的生态系统服务来测算其具体的服务价值, 即食物供、碳储量及生境质量, 这三项生态系统服务易受到城市化影响, 同时能够代表供给、支持和调节三种生态系统服务类型, 具体计算公式参照文献^[35]。

(1) 食物供给

不同土地利用类型生产食物的能力不同, 区域食物供给能力可以基于土地类型相应的食物能量转化来衡量^[35]。其中, 草地与奶类及食草动物肉类对应；耕地与粮食、油料、糖料、水果对应；水域与淡水产品对应(由于研究区未涉及到海域, 因此海水产品在本文中未被考虑)。通过ArcGIS 10.2重分类赋值到相应土地利用栅格数据上, 为使计算方便, 本文中食物供给值的单位为kJ/hm²。

(2) 碳储量

InVEST模型的原理是使用每种土地利用类型的碳密度及土地利用图来评估每个单元中的碳储存量^[35], 本文利用InVEST模型进行固碳量估算, 并考虑三个碳库, 分别为地上生物碳库、地下生物碳库和土壤有机碳库。通过ArcGIS重分类赋值到相应土地利用栅格数据上, 碳储量值的单位为t/hm²。

(3) 生境质量

生境质量是通过InVEST模型中Habitat Quality模块进行测算的^[35]。选取胁迫因子并对其权重进行赋权, 其权重大小体现了对生境类型的干扰强度, 干扰强度随着距离的增加而减小, 因此还需要设置每个胁迫因子的最大作用距离。同时, 不同生境类型对胁迫因子的响应都有所差异, 不同敏感性程度, 见表 3和表 4。

1.3.4 面板数据模型

生态系统服务受到多项因素的综合影响, 而这些因素与城乡建设用地关联有直接或间接的联系^[35]。因此, 本文选用面板数据模型来分析城乡建设用地关联(CLC)与生态系统服务价值之间的非线性关系^[40]。通过选取归一化差异植被指数(NDVI)、夜间灯光(NL)、人口密度(POPD)、GDP密度(GDPD)和交通(Traffic)作为控制变量。夜间灯光、人口密度、GDP密度和交通能够表现人类活动轨迹, 进而可以反映人类活动对生态系统服务变化的影响因素；NDVI代表自然条件, 可以通过影响人类活动强度和空气净化能力而间接影响城乡建设用地关联对生态系统服务的作用强度。面板数据模型为：

(4)

构建面板数据模型时, 为了避免非平稳性和异方差的影响, 所有变量都进行了自然对数变换；此外, 通过面板根单位检验和面板协整检验验证了面板数据模型的有效性。

2 城乡建设用地关联时空分析 2.1 城乡建设用地规模的时空变化

从2000—2018年城乡建设用地规模变化来看, 洞庭湖地区城镇建设用地规模呈现上升趋势, 农村居民点用地规模呈现下降趋势, 如表 5所示。其中城镇建设用地规模由140.43 km²增长至264.76 km², 增幅达到3.69%；农村居民点用地规模由1735.43 km²下降至1690.09 km², 降幅达到0.11%。

2.2 城乡建设用地关联强度的时空变化

洞庭湖地区城乡建设用地关联强度呈不断增强趋势。运用ArcGIS 10.2 Natural Breaks(Jenks)方法将关联强度值划分为5种类型, 即弱关联(0—1)、低关联(1—5)、中关联(5—10)、高关联(10—20)和强关联(>20)。从图 2可以看出, 研究期间城乡建设用地的高和强关联的网格数量不断增加, 其中, 强关联和高关联数量分别增加至6个和21个, 说明洞庭湖地区城镇建设用地与农村居民点用地的相互作用关系日益增强。从空间上看, 高关联和强关联的网格主要位于武陵区、鼎城区、资阳区、赫山区、岳阳楼区等市区范围, 呈现以市区为核心向周边扩张的趋势；弱关联和低关联的网格主要分布在距离市区较远的区域。主要是由于随着城镇化和工业化的发展, 市区周边的区域受经济社会辐射影响相对较大, 从而驱动周边区域经济发展和农村居民点用地规模的增长。

3 生态系统服务演变时空分析 3.1 生态系统服务演变时间变化

从各项生态系统服务价值的时序变化情况(表 6)可以看出, 生境质量和碳储量呈现下降趋势, 食物供给则呈现上升趋势。单位面积食物供给价值由599.47 kJ/hm²增加至876.28 kJ/hm²。单位面积碳储量价值由2067.96 t/hm²下降至2051.55 t/hm²。单位面积生境质量价值由0.648降低至0.635。

3.2 生态系统服务演变空间变化

从图 3可以看出, 各项生态系统服务价值存在明显的空间异质性。食物供给服务高值区主要分布在洞庭湖、沅水、澧水等水域周围的广阔沉积平原地区, 低值区主要集中在洞庭湖水域及研究区域边缘的山地丘陵地区。碳储量高值区主要分布于河湖近岸的湿地和研究区西南部及东部的山地丘陵地带, 低值区主要分布在河湖水域和城市建成区。生境质量较高的区域集中于洞庭湖水域及周边, 呈块状分散分布于东部、西南部地区, 较低的区域主要集中于城市建成区以及大片的山前平原农区。其中, 东部和西南部地区的碳储量和生境质量均出现明显下降趋势, 高碳储量区和高生境质量区明显减少。

4 城乡建设用地关联对生态系统服务的影响分析

如表 7所示, 所有变量在研究期间都是平稳的, 并相互协整。因此, 可以构建面板模型探讨城乡建设用地关联与生态系统服务之间的关系。

由表 8可知, 面板模型的校正决定参数Adjusted R²大于0.79, Prob (F-statistic)值小于0.001, 说明模型具有较高的拟合优度, 能够准确解释城乡建设用地关联与生态系统服务之间的关系。

城乡建设用地关联指数对生态系统服务有负向影响, 见表 8。城乡建设用地规模的增加导致了关联指数的变化, 从而造成对生态用地和农业用地的侵占, 直接导致生态系统服务价值下降。控制变量中, 人口密度和GDP密度对生态系统服务价值也呈现负向影响。城市扩张往往伴随着大量人口从农村向城市迁移和产业用地向城区外围布局。人口密度和GDP密度增加会导致一系列生态环境问题(如空气污染、土壤侵蚀、绿地退化), 这无疑会降低生态系统服务价值。夜间灯光和道路网络对生态系统服务有显著的负面影响。夜间灯光在一定程度上反映了人类活动的强度, 人类活动越多, 建设用地需求增加, 生态系统服务供给减少, 不利于生态系统的稳定发展。道路网络影响城乡建设用地关联强度, 改变区域景观格局, 导致区域生态风险增加, 对区域生态系统服务产生负向影响。NDVI与生态系统服务之间存在显著的正相关关系。NDVI越多, 植被覆盖度越高。植被覆盖度的增加可以提高生态系统维持自身结构、功能和抵御外部干扰的能力。

5 结论与讨论

以洞庭湖地区为研究对象, 结合2000、2005、2010、2015和2018年遥感影像提取的城乡建设用地信息, 综合运用引力模型、InVEST模型、面板数据模型等, 分析了洞庭湖地区城乡建设用地关联和生态系统服务的时空特征, 并探讨了两者间的关系, 结果表明：

(1) 洞庭湖地区城镇建设用地规模呈上升趋势, 农村居民点用地规模呈下降趋势。其中城镇建设用地面积由140.43 km²增长至264.76 km², 农村居民点用地面积由1735.43 km²下降至1690.09 km²。城乡建设用地关联强度呈不断增强趋势, 强关联的网格数量不断增长, 空间上呈现以建成区为核心向周边扩展, 低关联网格分布于离建成区较远的区域。

(2) 洞庭湖地区生境质量和碳储量呈现下降趋势, 食物供给则呈现上升趋势。单位面积食物供给价值由599.47 kJ/hm²增加至876.28 kJ/hm², 碳储量价值由2067.96 t/hm²下降至2051.55 t/hm², 生境质量价值由0.648降低至0.635。空间尺度上, 各项生态系统服务价值均存在明显的空间异质性。食物供给服务高值区主要分布在水域周围的广阔沉积平原地区, 低值区主要集中在洞庭湖水域及研究区域边缘的山地丘陵地区；碳储量高值区主要分布于河湖近岸的湿地和研究区西南部及东部的山地丘陵地带, 低值区主要分布在河湖水域和城市建成区；生境质量较高的区域集中于洞庭湖水域及周边, 同时在东部、西南部地区呈块状分散分布, 较低的区域主要分布于城市建成区以及平原农区。

(3) 社会经济和自然因素以及城乡建设用地关联对生态系统服务存在显著影响。其中, 城乡建设用地关联对生态系统服务有负向影响, 道路网络、夜间灯光、人口密度和GDP密度对生态系统服务也呈现负向影响, NDVI与生态系统服务之间则存在显著的正相关关系。

国家空间治理现代化背景下, 优化城乡空间发展格局是国土空间规划核心任务之一。城乡建设用地作为城乡经济社会发展的空间载体, 其关联强度的增加有利于构建协调有序的国土空间发展格局。本文研究表明洞庭湖地区城乡建设用地关联强度呈不断增强趋势, 体现出城乡空间要素流动更加紧密。同时, 城乡建设用地关联强度是生态系统服务下降的直接影响因素, 人口和经济的增长是间接影响因素, 其影响滞后于建设用地扩张。因此, 洞庭湖地区应该构建统一的城乡建设用地市场, 提高城乡空间要素流动的效率, 实现高质量的城乡空间融合, 而不是无序的扩大建设用地规模。同时, 建议加快卫星城市和小城镇的建设^[40], 进一步提升三个市区的辐射效应, 构建“刚性和弹性”边界管控体系, 减少对耕地或生态用地的占用, 促进整个洞庭湖地区国土空间高质量发展。同时, 洞庭湖地区生态系统服务价值的空间格局与该区域城市化建成区、农产品主产区和生态功能区的空间分布具有相似性^[41], 其中, 城市化地区城镇建设用地比例高, 其生态服务功能较弱, 形成了生态服务价值的低值集聚区。农产品主产区耕地面积广, 承担着主要的粮食生产功能, 形成了食物供给的高值集聚区。生态功能区以自然生态系统为主, 主要承担生态保护、生物多样性维持等功能, 城乡建设用地比例极小, 进而形成了生态服务价值的高值区。

综上所述, 本文对洞庭湖地区城乡建设用地关联和生态系统服务之间的关系进行了分析。但研究仍存在以下不足：首先, 本研究仅选取了三类生态系统服务来评估生态系统服务价值的变化, 后续的研究中, 需要选取更多类别的生态系统服务来获得更精确的生态系统服务价值。其次, 虽然面板数据模型已被证明可以有效地评价解释变量与自变量之间的关系, 但难以量化城乡建设用地关联与生态系统服务关系的空间依赖性和异质性。因此, 在接下来的研究中, 还需要利用空间回归模型(即空间滞后模型、地理加权回归模型和空间Durbin模型)进一步分析二者之间的关系。最后, 本研究仅在1 km尺度上探讨城乡建设用地关联与生态系统服务的关系, 未来的研究需要从多尺度视角来量化和比较城乡建设用地关联对生态系统服务的影响。