生态系统服务指人类从生态系统结构、过程和功能直接或间接得到的产品或惠益^[1]。以货币形式评估的生态系统服务价值(Ecological services value, ESV)是生态功能区划、环境保护及生态补偿的重要指标^[2-3]。土地利用/覆被变化(Land use and land use change, LUCC)是全球变化的重要驱动力, 通过影响生态系统格局与过程, 改变生态系统服务的提供能力, 最终导致生态系统服务价值发生变化^[4-6]。量化土地利用变化对区域生态系统服务的影响有助于优化配置土地资源和宏观调控、全面认识生态环境变化趋势并促进区域可持续发展^[7-8]。

国内外相关研究已从最初基于土地利用数量描述生态系统服务价值^[9-10]发展到应用多种指标、多种模型定量评价土地利用强度变化对生态系统服务的影响^{[7, 11]}, 并且随着精细尺度的研究需求不断增强, 研究者的关注范围也从全国^[12]、省级^[5]或区域^[6]尺度逐渐转向较小的空间尺度^{[4, 13-14]}。Chen等^[4]研究了1995—2015年全国县域生态系统服务强度与土地利用变化之间的联系；郭椿阳等^[13]基于格网探究了1990—2015年伏牛山区土地利用变化对生态服务价值的影响；雷军成等^[14]基于Landsat TM/OLI影像及当量因子法评估了1995—2016年江西省寻乌县生态系统服务价值的变化。县域是我国最基本的行政单元^[15], 近年来针对县域土地利用变化及生态系统服务价值方面的研究不断增多^[14-17], 但仍存在不足: (1)往往仅从数量上描述时间序列下土地利用变化造成的生态系统服务价值变化, 缺乏对生态系统服务价值损益进行空间分析；(2)仅对一级地类进行分析, 缺乏对同一地类内部生态系统服务质量差异进行描述；(3)多集中于生态系统服务价值现状评估, 不能预测未来的价值损益动态, 导致人们对未来县域精细尺度下生态系统服务价值损益的变化规律仍不清楚。

目前预测未来土地利用变化的方法主要有三种, 分别是经验统计模型、元胞自动机和系统动力学模型^{[16, 18-19]}, 其中CLUE-S模型(Conversion of land use and its effects at small regional extent)是发展较为成熟且应用最广泛的经验模型, 它基于土地利用与驱动要素的关系、对各种地类竞争及空间分配进行系统模拟, 适用于中小尺度研究, 而且具有较高的解释力与可信度^{[16, 20-21]}。该模型对县域土地利用变化的预测是量化未来生态系统服务价值损益时空动态的基础。

作为国家中心城市, 重庆市在城镇化推进过程中面临着日益严重的环境保护压力, 其中巴南区是重庆市主城九区中面积最大、城镇与农村发展相互作用最剧烈的区域^[22], 选择该区作为研究对象对了解县域尺度土地利用变化下生态系统服务价值损益的时空动态具有典型意义。基于2000年以来连续5期遥感影像及CLUE-S模型所预测出的未来10年巴南区8种土地利用类型的变化, 评估2000—2030年巴南区生态系统服务价值损益的时空动态、构建影响了强度模型度量土地利用变化对生态系统服务价值的影响, 旨在评估精细尺度下生态系统服务损益时空动态、提供科学合理的县域土地资源优化配置方案、推动城乡区域协调发展, 统筹中、小尺度国土空间规划并促进社会经济及生态的可持续发展。

1 研究区域与研究方法 1.1 研究区域

巴南位于重庆市西南部(106°26′2″—106°59′53″E、29°7′44″—29°45′43″N), 东接重庆市涪陵区和南川区, 南连綦江区, 西邻江津区、九龙坡区、大渡口区, 北傍南岸区、江北区、渝北区和长寿区(图 1)。巴南区东西横距46 km, 南北纵距70 km, 总面积1.83×10⁵ hm²。巴南区地貌复杂多样, 以丘陵和低山为主, 平地极少；地形起伏较大, 海拔范围为154—1133 m；年平均气温为18.7 ℃, 年均降水量为1000—1200 mm, 日照1100—1300 h；无霜期大于300 d, 属于亚热带湿润气候。全区生态环境良好, 森林覆盖率达47.6%, 拥有树种227种, 动物资源近百种, 并具有丰富的水资源与矿产资源。

1.2 数据来源

本文以2000—2020年连续5期土地利用/覆被数据为主要数据源。该数据来源于Landsat TM/ETM+/OLI遥感影像(http://glovis.usgs.gov/), 空间分辨率为30 m, 利用Erdas 9.2软件进行波段组合、几何校正、辐射校正、裁剪等预处理, 采用人机交互式目视解译。根据IGBP、FAO及马里兰大学的全球土地利用分类体系, 划分为水域、建设用地、未利用地、耕地和林地5种土地利用类型；同时结合巴南区土地利用/覆被特点和研究需求, 又将林地划分出针叶林、针阔混交林、灌木林和常绿阔叶林这4种二级地类。基于森林资源调查、固定样地等实测数据, 应用随机样点进行精度验证, 所有影像的总体分类精度均高于85%, 能够满足研究需求。海拔及地形数据来自地理空间数据云的SRTMDEM数据(http://www.gscloud.cn), 空间分辨率为30 m。该区城镇居民点、道路数据来源于国家基础地理信息中心(http://www.ngcc.cn/), 将这些作为预测未来土地利用格局变化的驱动要素。巴南区主要农作物类型、总产量、播种面积及价格等数据来源于《重庆市统计年鉴》^[23], 并结合2000—2015年《全国农产品成本收益资料汇编》^[24]估算生态系统服务价值当量。

1.3 CLUE-S模型应用 1.3.1 模型介绍

CLUE-S由荷兰瓦赫宁根大学研究团队开发, 基于每种土地利用类型与多个驱动因素的经验关系, 对研究区内所有地类进行明确的空间分配, 适宜于精细尺度土地利用变化的研究^[20-21], 主要步骤如下。

① 基于2000—2020年巴南区土地利用历史变化趋势, 应用马尔科夫链模型估算2025—2030年各类土地利用类型的面积。马尔科夫链模型是基于时间序列的随机过程, 广泛应用于地学数值模拟^[18]。根据巴南区现有土地利用类型集合X(t)预测未来各地类集合状态X(t+n), 基本形式为:

(1)

式中n为时间步长；P_ij为土地利用转移概率矩阵。

② 设置土地利用转化规则及转化弹性。根据巴南区现有土地利用变化规律设置土地利用转化规则, 将历史时段时内不同土地类型间发生过转化情况赋值为1, 将不同地类间未发生过转化情况赋值为0^[16]。弹性系数描述了不同地类转化的稳定性, 参考重庆市主城区已有的研究成果^[25]并经过多次模拟调试, 确定巴南区水域、建设用地、未利用地、耕地、林地(含所有二级地类)的弹性系数分别为0.9、0.7、0.5、0.5、0.6。

③ 选择地形因子(海拔、坡度、坡向)和区位因子(距道路距离、距城镇距离)作为驱动力对各地类进行Logistic回归(表 1), 得到各地类在空间某点的分布概率(P_{i, u})。

④ 根据分布概率、弹性系数以及迭代系数得到的总概率对各地类进行空间分配, 公式如下。

(2)

式中, TPROP_{i, u}为栅格i是土地利用类型u的总概率；P_{i, u}为Logistic回归求得的空间分布概率；ELAS_u为土地利用类型u的弹性系数；ITER_u为迭代系数。

巴南区8种土地利用类型与驱动因素的逐步回归结果表明所有土地利用类型均与坡度的关系最密切, 同时所有土地利用类型回归模型的ROC(Relative operating characteristics)曲线面积接近或大于0.7, 表明逐步回归模型能够较好地解释巴南区不同地类的发生概率。

1.3.2 马尔科夫链模拟检验及CLUE-S模拟结果检验

应用卡方检验对马尔科夫链所预测的未来各种土地利用面积结果进行验证。使用历史土地利用转移概率矩阵估算2020年巴南区各地类面积, 并与其真实土地利用格局进行比较(图 2), 应用Kappa系数检验各种地类空间分布格局预测的准确性。卡方检验结果和Kappa系数分别由公式(3)和公式(4)计算得到。

(3)

(4)

式中, O_i为观测值, E_i为模拟值, E_i为模拟期望；P_a为正确模拟的比例；P_c为随机情况期望模拟比例；P_p为理想情况下正确模拟比例。

根据公式(3)得到卡方值为0.0044, 小于χ_0.95(4)²=9.49, 即预测值与观测值差异不显著, 认为马尔科夫链预测的未来各种土地利用类型面积准确可靠。计算得到Kappa系数为0.8835, 即模拟精度为88.35%, 高于80%的检验标准。CLUE-S能够较准确地预测巴南区未来土地利用变化的空间格局, 结果合理可信。

1.4 生态系统服务价值评估

参考谢高地等^{[3, 26]}所提出中国陆地各种生态系统类型单位面积生态系统服务价值当量及生态经济学方法^[27](公式5—6), 构建了巴南区不同土地利用类型的各项生态系统服务价值当量表(表 2)。建设用地可提供的生态系统服务价值很低, 不参与计算^[28]。生态系统服务价值当量因子的经济价值为当年全国平均粮食市场价值的1/7^[26]。由《重庆市统计年鉴》^[23]和《全国农产品成本收益资料汇编》^[24]得到2000—2015年全国粮食平均出售价格范围为1098.80—2325.60元/t, 此期间全国与巴南区单位面积粮食产量范围分别为4.16—7.01 t/hm²、3.68—3.89 t/hm², 获得巴南区生物量修正因子为0.56—0.88。由于目前尚未能获得2020年粮食产量及价格数据, 因此在2020—2030年也以2015年的当量因子价值进行估算, 确保未来生态系统服务价值评估与历史时段结果具有可比性。同时为消除物价变动对生态系统服务价值评估的影响, 本文使用居民消费价格指数(Consumer price index, CPI)以2015年为基期, 对2000—2010年当量因子价值进行修正^[11], 得到2000—2030年巴南区生态系统服务价值的当量因子价值范围为818.53—1335.56元(图 3)。

1.5 土地利用变化对生态系统服务价值的影响

使用单一动态度和综合动态度分别描述各种土地利用类型及研究区整体土地利用变化的速率趋势及剧烈程度^[29]。基于土地利用转移矩阵核算得到的生态系统服务价值损益矩阵可揭示服务价值的数量变化特征、增减方向及来源^[7]；应用刘亚茹等^[7]所构建的影响程度指标来评价土地利用变化对服务价值的影响程度(公式7—10)。为从精细尺度分析2000—2030年巴南区生态系统服务价值损益的空间分布格局, 本文使用格网作为研究单元对巴南区生态系统服务价值损益做探索性空间分析, 了解损益的空间自相关特性。常见的格网尺度有500 m×500 m、1 km×1 km、3 km×3 km、5 km×5 km及10 km×10 km等^{[13, 16]}。考虑到本文研究尺度精细以及最小可塑性单元对结果的影响, 故而增加120 m×120 m和250 m×250 m这两种尺度, 通过分析不同格网尺度下生态系统服务价值损益时空动态及差异, 发现以250 m×250 m网格尺度的分析结果能够充分体现巴南区生态系统服务价值空间自相关特征, 是现有研究尺度、数据量和空间斑块情况下最合适的方案。因此本文将巴南区按250 m×250 m大小进行鱼网划分。按照自然断点法将生态系统服务价值损益从高到低分为5级, 使用Geoda软件进行价值损益的全局自相关分析分别计算全局、局部Moran′s I指数及局部G_i指数。

(5)

(6)

(7)

(8)

(9)

(10)

(11)

(12)

(13)

(14)

式中, V_{i, k}为生态系统服务当量价值(元)；P为粮食平均价格(元/kg)；Q为粮食平均产量(kg/hm²)；D_{i, k}为第i类生态系统第k项服务功能的当量因子；B为生态系统服务当量因子调整系数, 巴南区为0.671；ESV为生态系统服务价值(元)；A_i为土地利用类型i面积(hm²)；LP为生态系统服务价值损益(元)；VC_{i, k}、VC_{i, k′}分别为第i种土地利用类型的第k种生态系统服务在土地利用转变前和转变后的生态系统服务价值当量价值(元/a)；DY_i为单一土地利用动态度(%)；L_{i, t1}、L_{i, t2}分别为土地利用类型i在t₁年、t₂年的面积(hm²)；DY_t为综合土地利用动态度(%)；ΔL_i-j是t₁到t₂年期间第i类土地利用类型转化为非i类(j类)土地利用类型面积(hm²)；IM_i为土地利用变化对生态系统服务价值影响程度(元)；LP_{i, j}为第i类土地利用类型转化为第j类土地利用类型后生态系统服务价值损益(元)；I、I_i分别为全局和局部Moran′s I；G_i为局部G_i统计量；x_i为网格i在2000—2030年的生态系统服务价值损益值(元/hm²)；x是所有网格内生态系统服务价值损益平均值(元/hm²)；w_ij为网格i与网格j之间的空间权重, n为网格个数；S为网格i的生态系统服务价值损益方差。

2 结果分析 2.1 巴南区土地利用变化

2000—2020年巴南区耕地面积由1.39×10⁵ hm²减少了38.77%(图 4)。减少的耕地76.32%转化为林地、17.97%转化为建设用地；转化为水域和未利用地的比例最低, 分别为4.82%和0.88%。由于耕地的转化, 林地面积增长最多(4.1×10⁴ hm²), 未利用地增加最少(472.0 hm²)。CLUE-S模型预测上述土地利用变化趋势将在未来10年内延续。2020—2030年耕地继续减少1.85×10⁵ hm², 其中63.01%转化为林地, 29.62%转化为建设用地。2000—2030年耕地占巴南区面积的比例由75.99%降低到36.41%, 而林地占比增长了28.86%并将在2025年后取代耕地成为最主要的土地利用类型, 其中针叶林、灌木林、针阔混交林和常绿阔叶林的面积占巴南区面积比例增长依次降低, 分别增长了13.42%、7.15%、7.07%和1.22%(图 4)；建设用地和水域的比重分别增长8.29%和2.25%；未利用地占比最小, 其所占比重仅增长0.18%。

2000—2010年是巴南区土地利用变化最为剧烈的时期, 土地利用综合动态度由0.67%增加到1.08%(图 5), 主要来源于建设用地和林地的迅速增长, 其中建设用地单一动态度最高, 由12.16%增加到29.29%(图 5), 林地单一动态度也较高(4.94%—6.40%), 说明建设用地和林地是转化最显著的两种土地利用类型。此后巴南区土地利用变化速率放缓。2020年以后所有土地利用的单一动态度均趋于平缓, 综合动态度在2025—2030年最低(0.48%)。2000—2030巴南区建设用地、水域、林地、未利用地和耕地的单一动态度依次降低, 分别为41.58%、5.26%、4.55%、0.83%和-1.74%；耕地的单一动态度始终小于0, 即耕地持续转化为其他土地利用类型。

2.2 巴南区生态系统服务价值变化

2000—2020年巴南区生态系统服务价值由11.69亿元增长到27.64亿元(表 3), 此后仍逐渐增长, 到2030年达到30.19亿元。耕地生态系统服务价值在2010年达到最大值5.71亿元, 此后逐渐降低至3.38亿元。2000—2030年林地生态系统服务价值共增长16.19亿元, 其中针叶林和针阔混交林是巴南区生态系统服务价值增长最多的两种地类, 分别增加7.61亿元和5.06亿元；水体生态系统服务价值增长3.40亿元, 未利用地生态系统服务价值变化极小。2000—2030年耕地及林地的面积变化对生态系统服务价值影响最明显(图 5), 平均影响程度分别为32.24%和30.18%；水域面积增长对服务价值的影响逐渐增强, 而未利用地则基本无影响, 二者影响程度分别为5.60%和1.29%；建设用地对服务价值的影响程度为-2.39%, 表明其面积增加将导致生态系统服务价值的降低。2000—2030年除食物生产的服务价值将减少外(减少量为0.08亿元), 其他服务价值均增加(表 4)。气候调节价值增量最高, 为4.82亿元, 而维持养分循环价值增量最低, 为0.11亿元。耕地面积减少引起了食物生产价值持续下降, 而林地、水域面积的增长则是其他生态系统服务价值增长的主要原因。

2.3 生态系统服务价值损益时空动态

2000—2030年巴南区生态系统服务价值损失0.58亿元、总增益19.09亿元, 净增益为18.50亿元。耕地、林地转化为建设用地是生态系统服务价值损失的主要来源, 分别占全区损失的80.64%和8.42%(表 5)；而耕地转化为林地是服务价值增益的最主要方式, 贡献全区总增益的57.84%。土地利用结构变化和当量因子价值升高共同导致了生态系统服务价值增长, 二者分别贡献了净增益的70.96%和29.04%。

2000年巴南区耕地占绝对优势, 使全区生态系统服务价值呈现较为均匀的空间分布格局, 仅北部水系及海拔较高的桃子荡山的生态系统服务价值较高。随着时间推移, 城镇化发展的压力则使地势平坦的中西部地区耕地转化为建设用地, 使生态系统服务低值区由最初零星分布于西部边缘逐渐向中部及北部扩张。同时该区中南部、东南部海拔较高的山区耕地大面积转化为以针叶林为主的林地, 使生态系统服务价值明显增加, 使全区生态系统服务价值分布格局差异逐渐显现(图 6)。探索性空间分析结果显示(表 6), 巴南区生态系统服务价值损益的全局Moran′s I指数及P值分别为0.6077、0.001, 说明该区生态系统服务价值损益表现为显著的正向空间自相关。若将巴南区由东北到西南划一条对角线, 其西侧海拔较低、地势平坦而东侧海拔较高。生态系统服务价值冷点集中在对角线西侧海拔较高的李家沱、花溪街道、南泉街道等地区, 表现为“低值-低值”集聚特征和冷点(图 7)；而生态系统服务价值增益热点主要分布于对角线东侧海拔较高的石龙镇、石滩镇、接龙镇、跳石镇南部及安澜镇西侧, 呈现“高值-高值”集聚特征。

3 讨论

本文基于土地利用变化评估并预测了县域生态系统服务价值损益时空动态, 为未来区域土地资源的合理规划及有效利用提供了科学依据, 对保护生态环境和维持区域社会经济发展具有现实的指导意义。过去20年间巴南区最明显的变化是耕地持续转化为其他土地利用类型, 基于该趋势马尔科夫链模型预测出2030年巴南区耕地面积和林地面积分别为6.64×10⁴ hm²、9.13×10⁴ hm², 符合《巴南区及龙洲湾街道等22个镇街土地利用总体规划(2006—2020年)调整方案》对永久基本农田保护红线(4.66×10⁴ hm²)和生态保护(1.97×10⁴ hm²)的面积规定, 表明按现有趋势发展未来巴南区土地利用变化预测合理。

基于CLUE-S模型对精细尺度土地利用格局的准确模拟, 不仅能够获得未来生态系统服务价值的数量变化, 还可直观显示出不同土地利用类型及同一地类内部二级地类的生态系统服务价值损益的空间分异。2000—2030年巴南区土地利用结构变化导致的价值增长为13.13亿元, 占净增益的70.96%；而当量因子价值升高则贡献了净增益的29.04%。巴南区生态系统服务价值损益与地形密切相关。生态系统服务价值在西部低海拔地区随着城镇化加剧以及耕地转化为建设用地而损失逐渐加剧, 而在中南部及东南高海拔地区则随着耕地转化为林地而明显升高, 使价值损益在全局和局域尺度上均呈现显著的正向空间自相关特性。与我国同期大部分城市发展导致生态系统服务价值持续减少的情况不同^{[7, 11]}, 尽管巴南区建设用地逐渐扩张, 但林地面积的大幅增长使该区生态系统服务价值增益远高于其损失。这不仅反映出过去十几年间天然林资源保护工程及退耕还林工程的实施能够有效缓解城镇化与生态环境保护的冲突, 而且也说明现有土地利用发展模式能够兼顾经济发展及生态环境保护两方面的需求。

县域是介于省市及乡镇之间的中间尺度, 优化县域空间布局可有效推动大中城市和小城镇的协调发展。近期在重庆市“一区两群”发展格局规划中巴南区是主城都市区之一, 同时也是“国家城乡融合发展试验区”之一, 依据本文所得结果可从以下3方面优化该区未来空间布局, 以构建科学合理的城镇体系: (1)地势较平坦西部及北部地区应提高已有建设用地的利用强度、合理规划城市用地边界以避免减损水域及林地等生态系统服务价值高值区、促进产业结构调整和升级以增强中心城市综合承载能力；(2)南部海拔较高的地区继续进行退耕还林等生态保护工程、推行造林绿化、加强相邻区县间生态环境保护的联防联治、促进区域生境质量的全面提升；(3)由于耕地面积持续减少会引起食物生产服务的减弱, 因此需要严守耕地红线、使用先进技术提高粮食产量和保障粮食安全、搭建城乡产业协同发展平台, 推动城乡区域协调发展。

本文局限性主要包含以下3个方面: (1)CLUE-S难以反映政策对土地利用的直接影响, 所使用的Logistic逐步回归方程所拟合的土地利用类型与驱动因素的关系是单向作用过程, 不能反映土地利用类型变化的反馈机制与过程^[30]；(2)由于可获得数据及方法的限制, 本文未能预测未来生态系统服务当量因子经济价值的变化, 仍延用2015年的经济价值进行计算会造成未来生态系统服务价值低估或高估；(3)仅按照当前发展趋势对未来土地利用格局进行预测, 未能设置多种土地利用情景进行对比评估。今后需要应用更为准确的生态系统服务价值核算体系进行精细评估, 同时需依托各类生态系统定位观测研究站获取连续的自然驱动要素与社会经济数据, 设置多种发展模式进行模拟, 加强对微观土地利用变化过程与机理的研究、增强对生态系统服务流以及多种生态系统服务权衡与协同关系的认识。

4 结论

基于土地利用变化及CLUE-S模型可以直观评价并预测未来县域尺度生态系统服务价值损益时空动态。2000—2030年巴南区生态系统服务价值损失0.58亿元、总增益19.09亿元, 净增益18.50亿元；土地利用结构变化和当量因子价值升高分别贡献了净增益的70.96%和29.04%。耕地、林地转化为建设用地是导致巴南区生态系统服务价值损失的主要途径, 而耕地转化为林地则促进生态系统服务价值增益。巴南区生态系统服务价值损益与地形密切相关, 呈现显著的正向空间自相关特性。巴南区现有土地利用发展模式能够兼顾经济发展及生态环境保护两方面的需求, 但未来仍需实施生态工程、合理规划城市用地边界、严守耕地红线以实现区域长期可持续发展。