中国土地利用与生态服务价值空间动态研究——以地级及以上城市为例

危小建, 赵 莉, 程朋根, 谢亚娟, 王慧敏

(1.东华理工大学 测绘工程学院, 南昌 330013; 2.自然资源部 环鄱阳湖区域矿山环境监测与治理重点实验室, 南昌 330013; 3.仲恺农业工程学院 城乡建设学院, 广州 510225)

土地利用是人与自然进行紧密交流的重要纽带，土地利用/土地覆被变化影响着生态系统的结构和功能[1]，带来了生态系统服务价值的变化。生态系统服务(Ecosystem Services，ES)是指生态系统与生态过程直接或间接产生的，保证人类赖以生存的不可缺少的生命支撑产品和服务[2]，它是人类生存与现代文明发展的重要保障，同时也是科学技术无法取代的产物[3]。土地利用与生态系统相互作用、密不可分，分析土地利用与生态服务价值两者的相互关系，对合理规划土地资源，保证生态系统的协调发展具有重要意义。

近年来，国内外不少学者从多方面对土地利用与生态服务价值的关系开展了一系列研究[4-9]，已然成为生态学、地理学等领域的研究热点之一[10-11]。从研究驱动因素来看，茆长宝[12]、韩蕊[13]等基于土地利用数据和相关社会经济统计数据分别深入探究了促使江苏省和川东地区土地利用及生态服务价值产生变化的驱动因素；从研究方法和手段来看，张鹏岩等[14]采用热点分析、重心模型以及地理探测器模型研究了黄河流域生态服务价值的空间分异特征；姚小薇等[15]应用空间自相关理论分析了人口、空间、经济、生活4个不同层面城镇化水平对生态系统服务价值的影响。从研究尺度来看，目前开展的研究中主要是以单个省市县[16-18]、流域[19]或生态经济区[20]为研究对象，而针对中观尺度研究鲜少，以地级城市为研究单元也更为少见。同时，研究不同地级城市土地利用、社会经济发展与ESV的空间相关性，能够进一步揭示人类活动与社会经济发展对不同地区生态环境的影响，而当前从土地利用与社会经济两方面共同出发，对二者与生态服务价值的空间相关性方面也少有研究。因此，本文基于1990—2018年中国地级及以上城市土地利用数据和社会经济数据，分析土地利用与生态系统服务价值的动态变化情况，同时对生态服务价值进行敏感性分析，最后再通过研究生态服务价值与土地利用动态度、社会经济指标综合值的空间相关性，为优化土地利用结构、保证中国土地资源可持续利用以及促进生态文明社会协调发展提供一定的决策支持。

1 数据与方法

1.1 数据来源

本研究采用的中国地级及以上城市1990年、2000年、2010年、2018年的4期空间分辨率为1 km×1 km的土地利用数据由资源环境科学数据中心(http:∥www.resdc.cn)整理而得，土地利用类型主要通过耕地、林地、草地、水域、建设用地和未利用土地6种一级地类的形式呈现。相应年份的社会经济数据主要来源于中国统计年鉴和中国城市统计年鉴。本文研究对象是中国地级及以上城市单元，包括除香港、澳门和台湾地区以外的地级及以上城市和省直辖县级行政单元。西藏自治区无各个地级市GDP、城镇及农村人口指标等详细数据，本文以整个西藏地区的各项社会经济指标的平均值替代。另外，部分研究单元如镇江市、淮阴市及宿迁市等社会经济数据缺失，新疆维吾尔自治区无1990年城市统计年鉴等情况，使得部分数据难以获取，本文以相邻年份统计年鉴或同省周边相邻市县的数据插值替代。

1.2 研究方法

1.2.1 土地利用动态研究方法 土地利用动态度是分析土地利用结构变化动态的重要指标，主要包括单一土地利用动态度和综合土地利用动态度[14,21]。

(1) 单一土地利用动态度。单一土地利用动态度(K)是指在一定时间范围内某种土地利用类型的数量变化情况。其公式为：

(1)

式中:Ua，Ub分别为研究初期、末期某种土地利用类型的面积(km2)；T为研究时段。

(2) 综合土地利用动态度。综合土地利用动态度(LC)是指在某一时段所在研究区内，结合所有土地利用类型的变化情况而计算的综合土地利用变化的年变化率。其公式为：

(2)

式中：Ui为研究初期第i类土地利用类型的面积；ΔUi-j为研究时段第i类土地利用类型向其他非i类土地利用类型转化的面积的绝对值；n为土地利用类型的种类。

1.2.2 生态系统服务价值评估模型构建 Costanza等[22]在效用价值理论和均衡价值理论的基础上，首先对全球生态系统服务价值进行了评估，随后谢高地等[23-25]通过分析Costanza等[22]得出的全球生态系统服务价值，结合中国生态系统结构特点，进一步制定了中国生态系统单位面积的生态服务价值当量。本文借助谢高地等[24]的“中国陆地生态系统单位面积生态服务价值当量表”测算中国地级及以上城市的生态服务价值。其中，生态服务价值当量因子的经济价值约为全国平均粮食单产价值的1/7。参照前人研究[6,15,22]，建设用地的生态系统服务价值评估过程中需要考虑因素众多，测算较为复杂，本研究对其将不作估算。未利用地的生态服务价值估算时，其当量因子参数以草地与荒漠当量因子的平均值为基础进行估算[26]。因此，根据中国各省、自治区、直辖市不同粮食单产的实际情况计算当地单位面积粮食产量的经济价值，计算公式如下：

(3)

式中:Ea为研究区单位面积粮食产量的经济价值(元/hm2);P为2018年全国粮食平均价值(元/kg),本研究以国家发改委价格监测中心2018年的全国早中晚籼稻最低收购价格平均值与小麦最低收购价格的平均值作为2018年的粮食平均价格,取值为2.403元/kg;Q为根据1990年、2000年、2010年、2018年中国统计年鉴中31个省、自治区、直辖市粮食产量数据计算得到的1990—2018年各研究区域平均粮食单产(kg/hm2)(表1)。

表1 1990-2018年各地区粮食平均单产 kg/hm2

然而，随着社会不断发展，不同时间各研究区域社会发展程度有所不同，生态系统服务价值也会随之产生变化。为提高评价模型的精确性，本研究将以2018年为基准年，引入与社会发展相关的支付意愿系数(Dt)，构建适用于中国地级及以上城市的生态系统服务价值评估模型[27]，其公式如下：

(4)

式中:ESV为研究区生态系统服务总价值(亿元);Ai为第i类生态系统的面积(km2);Mij为第i类土地第j种生态系统服务功能所对应当量因子表数值;m为生态系统类型数;n为生态系统服务类功能项数;Ea表示单位面积粮食产量的经济价值(元/hm2)；Dt为引入恩格尔系数通过Logistic回归模型修正后的社会发展阶段系数,详见参考文献[27]。

1.2.3 生态系统敏感性指数 敏感性指数(Coeficient of Sensitive,CS)是指当生态服务价值系数变动1%时所引起的生态系统服务价值的变化情况[28]。其公式如下：

(5)

式中:ESVi和ESVj分别为调整1%前、调整1%后的生态服务价值;VCik和VCjk则分别为调整前后的生态服务价值系数。

1.2.4 空间自相关分析 空间自相关主要是用来研究某一变量在一定研究范围内是否成相关特征的空间分析方法，主要包括全局空间自相关(Moran′sI指数)和局部空间自相关(Local Moran′sI指数)[15]。

全局Moran′sI指数计算公式如下：

(6)

对于局部空间单元i，其Local Moran′sI指数计算公式为：

(7)

为了揭示多个要素的空间相关性，相关学者在Moran′sI指数的基础上进一步拓展了双变量全局自相关和局部自相关，其公式为：

(8)

2 结果与分析

2.1 土地利用变化分析

根据各年份土地利用数据，利用公式(1)和公式(2)分别计算1990—2000年、2000—2010年、2010—2018年以及1990—2018年的单一土地利用动态度和综合土地利用动态度。如图1—2和表2所示，中国地级及以上城市不同土地利用类型在空间上存在着明显的区域分异特征，林地、草地分布面积较大。从整个时间跨度来看，除草地呈减少的趋势外，其他土地利用类型均呈增加的趋势。其中，建设用地在研究区内变化速率最大，达2.50%，整个研究期共增加10.72万km2；水域和未利用地次之，整个研究期内土地利用动态度为0.29%，0.40%；耕地和林地变化速率则相对较小，整体上呈增加趋势，分别增加了1.50万km2，2.33万km2；同时，整个研究时段内综合土地利用动态度大于0.15%的区域主要集中于东北、西北和西部方向。

图1 1990-2018年中国地级及以上城市土地利用类型变化情况

注：基于标准地图服务系统下载的审图号GS(2016)1585号的标准地图制作，底图未做修改，下图同。

从单个时间段跨度来看，草地的土地利用动态度一直为负值，其减少幅度在2010—2018年最明显，为-1.17%；1990—2000年，耕地的土地利用动态度主要表现为增加的趋势，在2000—2010年、2010—2018年两个时段内均以较小的幅度减少，减少幅度不断增加；林地的土地利用动态度在1990—2000年为负值，在2000—2018年为正值，且在各个时间段始终小于0.1%，变化幅度较不明显；水域一直保持增加的趋势，尤其是2000—2010年期间变化最为强烈；未利用地在1990—2010年时段内以减少的趋势变化，变化幅度很小，在2010—2018年土地利用的变化幅度略微增大，呈增加趋势；建设用地面积则在研究时段内一直保持增加的趋势，变化幅度逐渐强烈，这主要是由于人口增加、城市扩张等原因，造成建设用地面积也随之逐渐增加；同时，综合土地利用动态度在2010—2018年变化最为明显，相对于其他类型土地而言，草地变化面积占总变化面积比重较大，且主要分布于西部地区，使得综合土地利用动态度变化明显区域主要集中于西部地区。

表2 研究区1990-2018年单一土地利用变化年动态度 %

2.2 生态系统服务价值变化分析

从时间尺度来看，研究区在1990年、2000年、2010年、2018年的生态服务价值ESV分别为：2.354×1013，2.253×1013，2.264×1013，2.197×1013元，整体呈减少的趋势。由表3可知，整个研究期内，生态服务总价值减少了1.57×1012元。1990—2010年期间，耕地的生态服务价值增加幅度逐渐减小，在2010—2018年则呈减少的趋势。林地和草地是我国最主要的两种土地利用类型，是构成生态系统服务总价值的重要组成部分，林地的生态服务价值在生态服务总价值中占比最大，在1990—2018年整个研究时段内，林地和草地两者的生态服务价值均呈减少特征，其中，林地减少2.8×1011元，草地减少了1.47×1012元。与生态服务总价值相比，由于水域占总面积比例小、未利用地产生的生态效用较小，使得两者产生的生态服务价值较小，变化幅度较为稳定。

从空间尺度来看(图3)，生态服务总价值超过1.5×1011元的区域主要分布于西藏、新疆和内蒙古等西北地区一带，这些地区海拔较高，草地分布较为广泛。生态服务价值在5×1010～1.5×1011元的区域主要集中于东北、东南和南方一带，这是由于所处地区气候和地理条件原因，使得林地面积分布范围广泛，产生生态服务价值较高。中东部则主要以耕地为主，该地带生态服务价值普遍低于5×1010元。

表3 1990-2018年研究区土地利用类型生态服务价值的变化 1010元

图3 1990-2018年研究区土地利用类型生态服务价值在空间上的变化

2.3 敏感性分析

通过公式(5)，将整个研究区域的各项生态服务价值系数分别变动1%，计算不同研究时期各土地利用类型的生态服务价值敏感性指数(表4)。整个研究期内，不同土地利用类型的敏感性指数一直处于稳定状态，且均小于1，说明生态服务价值系数在研究区内存在合理。另外，研究区内各土地利用类型敏感性指数从高到低依次为林地>草地>耕地>未利用地>水域。以上结果表明，林地对生态系统服务价值的贡献最大，未利用地和水域贡献最小。因此，为了保证生态系统服务价值的稳定协调，应加大对林地的保护。

表4 生态系统服务价值的敏感性指数

2.4 空间相关分析

2.4.1 生态服务价值变动与土地利用变化动态度的相关性分析 采用Geoda软件，构建空间权重矩阵，得到1990—2000年、2000—2010年、2010—2018年以及1990—2018年4个时间段内生态服务价值变化量与综合土地利用动态度的全局Moran′sI值，其结果为0.11，-0.08，-0.27，-0.18。结果表明，除1990—2000年，其余几个时间段内生态服务价值变化量与综合土地利用动态度均呈负相关。同时根据z检验结果，绘制1990—2018年整个研究期间内生态服务价值与综合土地利用变化动态度的双变量局部空间自相关LISA聚集图。由图4可知(高—高型表示生态服务价值与综合土地利用动态度均表现为较高值；高—低型表示生态服务价值较高，综合土地利用动态度表现为较低值；低—高型表示生态服务价值较低，综合土地利用动态度表现为较高值；低—低型表示生态服务价值与综合土地利用动态度均表现为较低值)，高—低型主要分布在中部地区，低—高型主要分布于西藏地区，低—低型主要分布于南部，高—高型则零星分布。

2.4.2 生态服务价值与社会经济的相关性分析 通过收集1990年、2000年、2010年、2018年的中国地级及以上城市总人口数、城镇人口数、城镇化率以及GDP这4个指标，利用归一化和熵值法求取社会经济指标综合值，再结合生态服务价值进行空间自相关分析，得到4个时期研究区生态服务价值与社会经济指标综合值的全局Moran′sI值，其结果分别为-0.14，-0.09，-0.10，-0.13。结果可知，生态服务价值与社会经济之间存在显著的空间负相关性，且两者在研究期内z检验值均<-1.96。根据生态服务价值与社会经济指标综合值的双变量空间自相关LISA聚集图(图5)，1990—2018年研究区生态服务价值与社会经济的局部空间特点总体上较为统一。其中，高—高型主要分布于中部和东部偏北地区，该区域经济较为发达，且林地所占面积较大，其生态价值较高；低—高型主要分布于中、东部地区，这部分地区经济发展较快，对生态资源占用现象严重；低—低型主要分布于西南地区，该区域山多地少、交通相对闭塞、地理条件较差，其经济和生态条件均处于落后阶段；高—低型则主要分布于昌都、噶尔县、拉萨和喀什市等西部地区，这是因为西部地区草地面积分布范围广泛，生态服务价值较高，但其经济发展较为落后。

图4 1990-2018年生态服务价值与土地利用变化动态度的

图5 1990-2018年生态服务价值与社会经济的双变量空间自相关LISA集聚图

3 讨论与结论

(1) 1990—2018年，草地在整个研究期内一直呈减少的趋势，共减少了38.74万km2；耕地在1990—2000年表现为增加的趋势，在2000—2018年呈减少的趋势，减少幅度略微增大，但整个研究期内仍保持增加趋势，净面积共增加了1.50万km2；林地在各个研究阶段变化幅度不明显；水域的土地利用动态度一直为正值，整个研究期内增加了2.06万km2；未利用地在1990—2010年以减少的趋势变化，2010—2018年有所增加，整体仍呈增加的趋势；建设用地在整个研究期内变化最为强烈，土地利用动态度始终为正值，增加了10.72万km2；综合土地利用动态度随着时间变化也有所增加，在2010—2018年变化尤为明显。

(2) 1990—2018年，中国地级及以上城市生态服务价值较高的区域主要集中于西北地区，总体上呈减少的趋势，28 a间共减少了1.57×1012元，而林地和草地产生的生态服务价值是构成研究区生态系统服务总价值的主要部分，其在研究期内分别减少了2.8×1011元、1.47×1012元。

(3) 1990—2018年，各土地利用类型的生态服务价值敏感性指数在研究区内均表现为<1，研究期内不同土地利用类型敏感性指数由高到低依次为林地、草地、耕地、未利用地以及水域。其中，林地敏感性指数最高，产生的生态服务价值最大。

(4) 1990—2018年，生态服务价值与综合土地利用动态度整体上呈负相关，各年份生态服务价值与社会经济指标综合值呈负相关，且双变量LISA图空间特点较为一致。低—高型主要分布于经济发展较快的中、东部地区；高—低型则主要分布于草地面积广泛经济发展较为落后的西部地区。这表明社会经济的发展、土地利用格局的变化一定程度上影响着生态系统结构和功能的转变。当地政府部门应结合实际情况，通过加大对生态系统的保护力度，采取有效措施保证土地合理规划，促进社会—经济—生态和谐发展。

与此同时，文中仍存在一些缺陷。本文仅根据谢高地等[23-25]制作的中国陆地生态系统单位面积生态服务价值表测算中国各地级城市生态服务价值，未根据各地区实际情况制定出适用于各地区的陆地生态系统单位面积生态服务价值表，导致测算有所偏差，在今后需要加强这一方面的研究。