生态系统服务定义为：由生态系统提供的, 人类从中收获的福祉和利益^[1]。通常使用经济价值来评价生态系统服务, 以确定我们从自然中获得的利益, 使不同区域具备可比性。在全球变暖的现状下, 碳减排已经变为一个全球性的紧要话题^[2]。两者因土地利用类型与陆地表面特征和功能而有很大区别。经济发展通过频繁改变土地利用方式间接影响陆地生态系统和人为的碳排放。土地利用变化对碳排放不只是简单的线性关系, 而是更为深远的影响。同时土地利用方式的改变也极大地影响了自然生态系统, 从而影响了它们所提供的服务。因此土地利用变化不仅影响着碳排放的空间分布、范围和强度, 而且还通过改变生物多样性、区域资源和生态系统类型的空间分布间接影响ESV^[3]。随着工业化的快速发展, 中国成为全球碳排放量最高的国家^[4]。特别是1950—2005年因土地利用变化引起的累积碳排放量为106×10⁸ t, 是同期世界土地利用变化碳排放量的12%^[5]。为了减缓碳排放和全球变暖, 中国制定了不同的环境政策并采取了有效措施^[6]。特别是, 2020年习近平总书记在第七十五届联合国大会一般性辩论上提出了碳达峰、碳中和的目标。2021年国家林草局为落实这一决定提出“十四五”期间将扎实开展林业和草原碳汇行动, 提升林草生态系统碳汇能力。因此, 在碳减排政策和改善生态环境的双重背景下, 研究两者间的空间关系模式, 对解决社会发展过程中既要保持经济增长, 又要兼顾生态环境这一难题有重大意义。

自20个世纪90年代开始, 一方面, 国内外学者围绕土地利用碳排放问题从国家^[7]、省域^[8—9]、城市^[10]、区域^[11—12]等不同尺度对碳排放与经济增长的脱钩效应分析^[13—15]、单一地类碳排放机理剖释^[16—17]、影响因素^{[7, 18—19]}、变化规律^{[10, 20]}、效率^[21]、地类变化与碳源/汇的关系^[11]、碳排放空间差异与关联^[22]等角度展开了大量的研究。另一方面, 在评估土地利用变化造成ESV波动方面, 学者们还进行了许多研究^[23—26], 如ESV和土地利用的时空变化特征^[27—28]、相互作用关系^{[25, 29—30]}等。综上可知, 以往都是从不同角度和尺度对两者单独开展了研究, 但是关于土地利用碳排放和ESV之间空间关系的研究较少^[31—32], 并且南四湖流域还未开展土地利用碳排放的相关研究。

南四湖流域位于鲁南经济圈, 是我国北方最大的淡水湖泊, 肩负重要粮食生产的重任, 同时在京杭大运河和国家战略性工程中发挥重要作用。未来如何推动流域内各县市在经济持续发展的前提下, 实现节能减排、保护生态环境是值得思考的问题。鉴于此, 本文以南四湖流域为例, 利用土地利用数据和统计资料, 对2000—2018年间的土地利用转移中碳排放量和ESV的时空演变开展研究, 并基于格网对碳排放强度和生态系统服务强度进行空间关联分析, 以期为流域提出碳减排政策, 建立低碳土地利用结构, 努力改善生态环境总体质量, 从而走向绿色低碳发展的道路提供支撑。

1 研究区概况

南四湖流域处于中国东部淮河流域的北部, 涉及山东、江苏两省5市, 面积约3.17×10⁴ km²(图 1)。地处鲁中南山区西侧的山前冲击洪积平原, 中、西部为平坦的平原, 东部为山地丘陵地形。海拔在-42—630 m之间, 绝大部分地区低于70 m。耕地分布广泛, 占全区面积的68.5%以上。自21世纪以来由于城镇化进程的推进, 流域国内生产总值由1174.91×10⁸元变为12078.12×10⁸元, 翻了9.28倍；建设用地面积从4635.15 km²增长到6035.81 km², 增长了30.22%。近年来, 该区受社会发展的直接影响, 土地利用方式剧烈变化, 能源消耗需求增加, 对当地的生态和碳排放系统带来一系列问题。同时由于其所处地理位置及其作用, 各级政府和相关部门都十分重视生态保护, 启动了《南四湖流域水污染综合整治三年行动方案(2021—2023年)》等一系列工作, 以防止环境恶化。

2 数据来源与研究方法 2.1 数据来源及预处理

本文使用地理国情监测云平台的2000—2015年的土地利用数据, 其分辨率为30 m。2018年土地利用数据来源于中国科学院资源环境科学数据中心, 是通过人工目视解译Landsat 8遥感影像获得的。依据中科院土地利用分类标准, 将研究区的地类分为6个一级类型, 依次为耕地、林地、草地、水体、建设用地和未利用地。能源消耗和社会经济数据分别来源于对应年份的《中国能源统计年鉴》、《统计年鉴》和《全国农产品成本收益资料汇编》等统计资料。利用网格法按照5 km×5 km大小的网格对研究区进行划分, 共得到1387个评价单元。由于数据可获取性, 本文中2000—2018年南四湖流域能源消耗数据, 根据研究区内各区县GDP与所在省份GDP的比值来进行推算。

2.2 研究方法 2.2.1 土地利用类型转移矩阵

土地利用转移矩阵以二维矩阵的形式反映19年间各用地类型的数量转移情况, 清楚掌握用地类型在2000年的流动方向及2018年的来源构成。其公式为：

(1)

式中：S_ij为土地面积；n为6；i、j分别为2000年与2018年的地类序号。

2.2.2 碳排放量的估算

研究区碳排放量计算：即各用地类型的面积和相应的碳排放系数相乘再求和。碳排放计算公式为：

(2)

式中：E_Κ表示直接碳排放量；e_i表示各用地类型的碳排放量；S_i和δ_i分别表示用地类型i的面积和碳排放系数, 参考相关文献^[33—34]和研究区的现实情况确定各用地类型的碳排放系数, 见表 1所示。

基于2000—2018年南四湖流域化石能源消耗量计算建设用地碳排量, 然后除以其面积得到建设用地碳排放强度。建设用地碳排放量的计算采取间接估算法, 结合《中国能源统计年鉴》和2006年IPCC《国家温室气体排放清单指南》(表 2)^[38], 通过各种能源的消耗量和碳排放系数估算得到。

建设用地碳排放计算公式：

(3)

式中, E_ŋ为建设用地碳排放量；e_i为能源i消耗产生的碳排放量；E_i为能源i的消耗量；μ_i为能源i消耗量转换为标准煤的系数；ε_i为能源i的碳排放系数。

土地利用碳排放强度由格网内各用地类型的面积和碳排放系数决定, 即碳排放强度越大说明单元格网的碳排放量越大。计算公式如下：

(4)

式中, C_RI为土地利用碳排放强度；S_i和P_i分别为格网内用地类型i的面积和碳排放系数；S为格网总面积。

2.2.3 生态系统服务价值

参照谢高地等^[39—40]的研究, 1个生态系统服务价值当量的价值量由每公顷农田年平均自然粮食产量的经济价值1/7^[40—41]来确定, 已广泛在流域^[39—40]等尺度开展深入研究。为避免各年份和不同省份粮食价格的差异对价值计算的影响, 本文仅以山东省2018年小麦、玉米和稻谷3种主要农作物的播种面积、单产和粮食价格为基础数据, 计算出南四湖流域ESV的单位当量为2085.53元, 其计算模型如下：

(5)

式中, E_a为南四湖流域ESV当量因子的经济量；n为3；p_i为农作物i在山东省的价格；q_i为农作物i单产；m_i为农作物i面积；M为农作物总面积。

由于我国疆域辽阔, 各地的生态系统各具特点, 为降低直接使用2015年谢高地完善的“单位面积生态系统服务价值当量表”^[39]产生的误差, 本文参考研究区已有研究^[42]和现实情况修正了各地类的生态系统服务价值系数(表 3), 其中不考虑建设用地的ESV。

生态系统服务价值：

(6)

(7)

生态系统服务价值强度：

(8)

式中, A_ij为网格j中地类i的面积；C_i为地类j的价值系数；S为网格面积；n为6种土地利用类型, m为1387。

2.2.4 双变量空间自相关性分析

我们使用双变量Moran′s I研究了碳排放强度与ESV强度之间的空间聚类和离散。本文使用了两种类型的双变量Moran′s I, 即全局双变量Moran′s I和局部双变量Moran′s I。全局双变量Moran′s I研究了整个区域内碳排放强度与ESV强度之间是否存在以及空间相关性程度, 而局部双变量Moran′s I则在具体的空间格网上显示空间相关性。计算公式如下：

(9)

(10)

式中, I和I_i分别为碳排放强度与ESV强度的全局和局部双变量Moran′s I。n是空间网格单位的总数。w_ij一个n×n空间权矩阵, x_i和x_j为单元网格i和j的属性值, x是属性值的均值, S²为方差。双变量LISA集聚图可视化局部空间相关性。它说明了特定位置的碳排放强度值与相邻位置的ESV强度值在一定显著性水平上的关系。双变量LISA生成的4个象限代表了4种类型的局部空间自相关。

3 结果与分析 3.1 土地利用变化分析

如表 4所示, 2000—2018年南四湖流域各土地利用类型间均进行了复杂的转移, 发生转移的总面积为4468.86 km², 是流域总面积的14.10%。耕地对转移面积的贡献最大, 为2579.12 km², 其中有2084.55 km²变为建设用地, 占比80.82%, 15.72%的转出量转为水体, 虽占比较低, 但转出面积大, 为405.40 km²；其次是建设用地, 转出量为854.38 km², 95.87%的土地转为耕地；其它大部分地类的转出量转为耕地。出现此情况的原因是流域经济快速发展促使城镇外围的耕地被大量占用, 但为保证粮食安全, 故改变其它生态用地的使用方式, 从而确保耕地总量, 之后又由于实施“两退三还”(退耕、退池、还林、还湿、还湖)工程和居民点整治等原因导致土地转换复杂且多样。建设用地是转入最多的土地类型, 转入2255.05 km², 其次是转入面积为1581.84 km²的耕地。其他地类均发生转入转出, 情况各异。

3.2 土地利用碳排放分析 3.2.1 南四湖流域2000—2018年土地利用碳排放测算与分析

研究区土地利用碳排放量的改变受地类和土地利用程度的影响。本节主要以各地类的碳排放量、碳源/汇和净碳排放量等指标测算南四湖流域土地利用碳排放, 分析2000—2018年间各指标的动态变化特征。利用上述计算公式和相关收集数据计算南四湖流域在研究期间的相应碳排放指标, 结果见表 5。

通过计算可知, 南四湖流域的净碳排放量表现为碳排放(表 5)。2000—2018年碳排放能力显现稳定增长趋势, 由1262.51×10⁴ t增长至6438.41×10⁴ t, 年均增长率为22.78%, 变化趋势与碳源量相似。其中, 建设用地碳排放量增长迅速, 从2000年的1176.36×10⁴ t增长至2018年的6355.35×10⁴ t, 但增长量具有“先快速, 后放缓”的特点, 这是由于流域内城镇化水平得到大幅度提高, 伴随着大量的能源消耗和建设用地增加, 之后由于能源结构的优化和节能减排的实施；耕地的碳排放量由于本身面积的减少呈现缓慢下降的趋势, 由2000年的96.02×10⁴ t下降至2018年的91.81×10⁴ t, 年均下降率为0.24%。碳汇量在19年间呈现先增加后减小并趋于稳定的态势。其中林地和草地的碳吸收量在2005年后出现下降, 水体基本保持不变, 未利用地波动较为明显, 呈先减少后增加的特点。

从碳源/汇构成上分析, 碳源方面, 建设用地的碳排放是最主要碳源, 占总碳源量的92.45%—98.58%, 并呈逐年增长态势。同时耕地作为碳源之一, 碳排放贡献率低, 但面积占比高。地均建设用地碳排放量表示各年建设用地的碳排放强度以及区域的能源消费结构, 自2000年以来表现为快速增长趋势, 各年的地均建设用地碳排放量依次为25.38、62.16、80.28、96.02、105.29 t/hm², 其中2018是2000年的4.15倍, 因此建设用地在碳减排方面具有较大潜力；碳汇方面, 林地和水体对碳吸收的贡献率相当, 两者之和占总碳汇的96.11%以上, 说明林地和水体具有重要的碳汇能力, 而在流域内两者占总面积的5.89%—8%, 凸显了其重要的碳吸收能力, 另外草地和未利用地的碳吸收量较小。由于碳源量与碳汇量相差悬殊, 因此该流域今后的工作重心是碳源的控制与碳汇的增加。

3.2.2 土地利用碳排放强度的时空特征

由于流域内地类结构和发展水平的差异, 因此为了具有可比性以及直观地分析土地利用碳排放量的变化情况, 故本文利用碳排放强度指标, 依次分为Ⅰ、Ⅱ、Ⅲ、Ⅳ和Ⅴ 5个区间(见图 2)。19年间南四湖流域碳排放强度时空变化明显, 其特点如下：流域碳排放强度低值多分布于东部, 整体上呈现增加的特点, 碳排放强度最大值从2000年的21.61 t/hm²增长到2018年的101.42 t/hm², 增长了4.69倍。落在Ⅰ区间的格网数逐年减少, 尤其是2000—2005年间减少了78.3%, 出现了向Ⅱ区间转移的现象, 而同时南四湖湖区始终属于Ⅰ区间；落在Ⅱ区间的格网呈现先增加后减少的特征, 2010年以来土地利用类型主要为林地和草地的东部和南四湖湖区的周边稳定的处于该区间；Ⅲ区间的格网数的增长率最高, 在2018年占总格网数的46.36%, 绝大部分与耕地的分布范围相吻合；Ⅳ和Ⅴ区间的格网从无到有并逐年增加, 呈团块状, 主要集中在人类活动活跃的建设用地。

3.3 生态系统服务价值的时空演变

2000—2018年南四湖流域生态系统服务总价值呈波动变化, 总体上呈增加趋势, 总价值704.07×10⁸元增加到743.37×10⁸元, 增加了5.58%(表 6)。从各地类的ESV来看, 耕地和水体是总ESV的贡献主体, 依次占24%和61%以上。19年间水体的ESV占比稳定增长, 由61.85%到67.79%；而耕地的ESV占比逐年下降, 由27.03%到24.48%；其它地类的ESV总体上呈下降态势。从土地利用构成中得到, 19年间建设用地的面积大幅度扩张, 林地、耕地、草地和未利用地均发生不同程度的减少, 而水体面积的增加, 反而使总ESV增加29.30×10⁸元, 说明水体的ESV强度最高, 具有极其重要的生态地位, 同时各地类间的相互转换是导致研究期间内总ESV波动变化的关键原因。

为了直观地分析流域ESV的空间分布特征, 同时避免流域边缘面积的影响, 本节利用ESV强度表示各格网, 并对其重分类依次分为Ⅰ、Ⅱ、Ⅲ、Ⅳ 4个区间(见图 3)。南四湖流域的ESV强度具有明显的空间差异, 西部地区的ESV强度整体高于东部地区。其中, 处于Ⅰ区间的ESV强度的分布范围广, 其格网占总格网的72.45%以上, 主要由于耕地的空间分布相对均一且占比高；处于Ⅳ区间的格网占总格网的4%左右, 土地利用类型为水体, 主要分布在南四湖湖区。与其他年份相比, 2018年处于Ⅰ区间的格网变少, 出现向Ⅱ区间转变的现象, 同时Ⅱ区间的分布范围主要位于流域的东部, 集中分布着草地和林地；Ⅳ区间的分布范围由原先的条带状出现断裂；Ⅲ区间的分布范围变化不明显。

3.4 碳排放强度和ESV强度的关系

为了深入地了解碳排放强度与ESV强度的相互交互作用, 首先利用SPSS软件得到2000—2018年两者的相关系数依次为-0.453、-0.413、-0.347、-0.335和-0.341；然后借助空间自相关揭示两者的空间关联规律(表 7), 各年份的Moran′s I均为负, P值小于0.001, 表明研究区内土地利用碳排放强度与ESV强度之间具有空间负相关性, 即随着碳排放强度的增强, ESV强度呈下降的态势, 并且Moran′s I具有逐年下降的特点, 表现出整体较为分散、小范围聚集的空间布局格局。

通过LISA集聚图可以看出(图 4), 2000—2018年间, 低高聚集区稳定地分布在南四湖湖区, 表示该区碳排放强度低值与ESV强度高值形成集聚现象；2000和2005年的高低聚集区集中分布在耕地面积占比高, 建设用地点状零散分布的西部和北部, 而在2005年之后出现了大面积向低低聚集区转变的现象, 高低聚集区与呈片状分布的主要城镇用地相吻合, 主要是由于2005—2010年间地类之间的变化剧烈, 总转移面积为3581.62 km², 其中有1687.82 km²的其他类型土地转化为建设用地, 说明在该时间段流域的城镇化水平得到快速发展, 能源消耗增加, 导致局部地区的碳排放强度发生大幅度的上升；虽然建设用地和耕地的价值系数低且都属于碳源, 但是前者的碳排放系数要远高于耕地, 从而使耕地面积占比高的格网与周围的低ESV强度形成集聚现象；低低聚集区在2000—2005年零散分布且变化不大, 之后分布范围扩大, 主要分布在西部和北部, 地类为具有面积优势的耕地；高高聚集区呈零星分布, 这是由于周围是水体, 属于碳汇。

4 讨论

以往的研究大多集中在独立的行政单元(如县、市等)的经济和社会等方面的碳排放效应分析^{[7—9, 18]}。本文将跨越行政区的流域作为一个研究整体, 基于格网探究因土地利用变化引起的碳排放和ESV以及两者的空间聚集规律。研究发现流域的建设用地面积和净碳排放量在这19年内呈增长的趋势, 与Chen等^[43]以成渝城市群为研究区、李玉玲等^[14]以陕西省为研究区的研究结果一致, 说明土地利用的改变是影响碳排放量和碳源/汇格局的关键原因, 不当的土地开发会影响碳平衡；同时与中原城市群^[31]的空间关系一样, 碳排放强度与ESV强度之间具有空间负相关性, 碳排放强度对ESV强度有空间溢出效应。结合研究结果, 通过在流域内实施能源消费结构的升级、提高利用效率、引进低碳生产技术、防范耕地被滥用和提高建设用地集约度, 可以实现碳排放的降低和有效保护生态系统。

尽管本研究展示了19年间碳排放和ESV的时空特征, 但仍存在一些潜在的局限性需进一步改进。首先, 本文忽略了空间异质性, 只是借鉴前人研究和资料得到各地类地碳排放/吸收系数^[44], 如耕地碳排放受种植作物种类、土壤的性质和耕作特点等多种因素的影响而有所差异。本文与国内绝大部分研究类似将耕地视为碳源, 而马涛^[45]的研究表明上海的农业活动表现为碳汇。各能源碳排放系数参考了2006年IPCC国家温室气体清单指南, 这可能并不完全适用于研究区。因此, 今后根据区域的实际情况科学合理地确定碳排放系数, 或者形成统一的测算方法；其次, 本文使用了高精度的土地利用数据, 但基于5 km×5 km网格对流域进行划分。由于生态和地理过程伴随着复杂的尺度效应^[46], 不同的尺度研究的侧重点不同。因此未来有必要从多尺度探究碳排放强度与ESV强度关系特征；最后, 由于受数据收集的限制, 本文通过GDP占比推算流域的能源消耗, 同时没有将工业生产、生活和废弃物排放加入建设用地的碳排放核算中, 所以计算结果有误差, 但不影响区域碳排放的时空分析。因此, 今后需要寻找全面且受数据限制小的计算方法。

5 结论

(1) 研究区具有面积优势的地类是耕地和建设用地, 分别占总面积的68.5%和14.5%以上。在19年内各地类之间进行了不同程度的转移, 其中耕地和建设用地是变化量最大的类型。产生该变化的原因与城镇化、人类生产生活和退耕还湖等相关政策有关。

(2) 19年内流域的总价值随地类之间的转化呈波动变化, 整体上是增加, 水体面积的增加是导致其的决定性原因。通过ESV强度图可知, 东部的ESV强度明显高于西部, 同时Ⅲ、Ⅳ区间分布在南四湖湖区, Ⅱ区间分布在海拔较高的东部, 这与地类、价值系数和地形有关。

(3) 南四湖流域的净碳排放量表现为碳排放, 且稳定增长。其中, 碳源量要远大于碳汇量, 地均建设用地碳排放量19年间增长了4.15倍, 这与经济社会发展有密切的关联。通过碳排放强度的时空演变分析可知, Ⅰ区间的格网转变为Ⅱ和Ⅲ区间的格网, 后者的分布范围迅速扩张, 同时Ⅳ和Ⅴ区间的格网从无到有且逐年增加, 与主要建制镇的位置相吻合。整体上呈现东部高、西部低的特点, 与ESV强度的空间分布有一定的关联。

(4) 对南四湖流域土地利用碳排放强度和ESV强度进行数量和空间统计分析, 得知两者具有负相关性, 且通过了P值检验。19年内两者的空间关联规律发生较大的变化。